API tools faq
paste
SpaceQuester

Untitled

SpaceQuester
Nov 10th, 2019
451
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 28.87 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.  
  11. #define Node_count 2000
  12.  
  13. #define Equations_per_node 8 // !!! Don't change !!!
  14. #define Equations_count Node_count * Equations_per_node
  15.  
  16. int MaxDeep = 3;
  17.  
  18. double f[Equations_count];
  19. double f_diff[Equations_count];
  20.  
  21. bool enable_I_syn_out = false;
  22.  
  23. double c_0 = 2; // uM
  24. double c_1 = 0.185;
  25. double v_1 = 6; // s^-1
  26. double v_2 = 0.11; // s^-1
  27. double v_3 = 2.2; // uM/s
  28. double v_4[Node_count]; // uM/s - Controling parameter // 0.5 //double v_4[Node_count]; // uM/s - Controling parameter //0.495
  29. double v_5 = 0.025; // uM/s
  30. double v_6 = 0.2; // uM/s
  31. double k_1 = 0.5; // s^-1
  32. double k_2 = 1; // uM
  33. double k_3 = 0.1;
  34. double k_4 = 1.1; // uM/s
  35. double a_2 = 0.14; // uM/s
  36. double d_1 = 0.13; // uM
  37. double d_2 = 1.049; // uM
  38. double d_3 = 0.9434; // uM
  39. double d_5 = 0.082; // uM
  40. double alpha = 0.8;
  41. double tau_IP3 = 7.143; // s
  42. double IP3_star = 0.16; // uM
  43. double d_Ca = 0.001; // 0.001
  44. double d_IP3 = 0.2; // 0.12
  45. double alpha_Glu = 2; // 2
  46. double g_astro = 3; // 3
  47.  
  48. // https://neuronaldynamics.epfl.ch/online/Ch2.S2.html
  49. double C_m = 1; // muF/cm^2
  50. double g_K = 35; // mS/cm^2
  51. double g_Na = 40; // mS/cm^2
  52. double g_L = 0.3; // mS/cm^2
  53. double E_K = -77; // mV
  54. double E_Na = 55; // mV
  55. double E_L = -65; // mV
  56.  
  57. double p_rewir;
  58. double p_inhib;
  59.  
  60. double I_app_min;
  61. double I_app_max;
  62.  
  63. double g_syn;// 0.18 // 0.04 // 0.2 // 1.6
  64. double k_syn = 0.2; // 0.2
  65. double E_syn[Node_count];
  66.  
  67. double alpha_G = 25; //s^-1
  68. double beta_G = 500; //s^-1
  69.  
  70. double I_app[Node_count];
  71.  
  72. double** A_A;
  73. double** B_A;
  74. double* C_A;
  75.  
  76. double** A_N;
  77. double** B_N;
  78. double* C_N;
  79.  
  80. FILE* fp_I_syn;
  81.  
  82. double tau[Node_count][Node_count];
  83.  
  84. #define tau_min 2 // ms
  85. #define tau_max 12 // ms
  86.  
  87. #define ms_to_step 40 // (0.001 / dt) !!! Don't forget !!!
  88.  
  89. #define Max_delay tau_max * ms_to_step
  90. double V_old_array[Node_count][Max_delay];
  91.  
  92. const double Freq = 500; // Hz
  93. const double Min_magintude = -0.13; // pA
  94. const double Max_magintude = 8.0; // pA
  95. const double Duration = 0.001; // sec
  96.  
  97. double Meander_start_from_zero[Node_count];
  98. double Meander_width[Node_count];
  99. double Meander_height[Node_count];
  100. double Meander_interval[Node_count];
  101. double last_meander_end[Node_count];
  102.  
  103. double I_stim(int i, double t)
  104. {
  105. if (t < Meander_start_from_zero[i])
  106. return 0;
  107.  
  108. t -= Meander_start_from_zero[i];
  109. t = fmod(t, Meander_width[i] + Meander_interval[i]);
  110.  
  111. return t < Meander_width[i] ? Meander_height[i] : 0;
  112. }
  113.  
  114. double Ca(int i)
  115. {
  116. return f[i * Equations_per_node];
  117. }
  118.  
  119. void SetCa(int i, double value)
  120. {
  121. f[i * Equations_per_node] = value;
  122. }
  123.  
  124. double IP3(int i)
  125. {
  126. return f[i * Equations_per_node + 1];
  127. }
  128.  
  129. void SetIP3(int i, double value)
  130. {
  131. f[i * Equations_per_node + 1] = value;
  132. }
  133.  
  134. double z(int i)
  135. {
  136. return f[i * Equations_per_node + 2];
  137. }
  138.  
  139. void Setz(int i, double value)
  140. {
  141. f[i * Equations_per_node + 2] = value;
  142. }
  143.  
  144. double G(int i)
  145. {
  146. return f[i * Equations_per_node + 3];
  147. }
  148.  
  149. void SetG(int i, double value)
  150. {
  151. f[i * Equations_per_node + 3] = value;
  152. }
  153.  
  154. double V(int i)
  155. {
  156. return f[i * Equations_per_node + 4];
  157. }
  158.  
  159. void SetV(int i, double value)
  160. {
  161. f[i * Equations_per_node + 4] = value;
  162. }
  163.  
  164. double m(int i)
  165. {
  166. return f[i * Equations_per_node + 5];
  167. }
  168.  
  169. void Setm(int i, double value)
  170. {
  171. f[i * Equations_per_node + 5] = value;
  172. }
  173.  
  174. double n(int i)
  175. {
  176. return f[i * Equations_per_node + 6];
  177. }
  178.  
  179. void Setn(int i, double value)
  180. {
  181. f[i * Equations_per_node + 6] = value;
  182. }
  183.  
  184. double h(int i)
  185. {
  186. return f[i * Equations_per_node + 7];
  187. }
  188.  
  189. void Seth(int i, double value)
  190. {
  191. f[i * Equations_per_node + 7] = value;
  192. }
  193.  
  194. double V_old(int i, int delay)
  195. {
  196. return V_old_array[i][Max_delay - 1 - delay];
  197. }
  198.  
  199. int RandomI(int min, int max)
  200. {
  201. return ((double)rand() / (RAND_MAX - 1)) * (max - min) + min;
  202. }
  203.  
  204. double RandomD(double min, double max)
  205. {
  206. return ((double)rand() / RAND_MAX) * (max - min) + min;
  207. }
  208.  
  209. double J_channel(double* f, int i)
  210. {
  211. 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);
  212. }
  213.  
  214. double J_PLC(double* f, int i)
  215. {
  216. return v_4[i] * (Ca(i) + (1 - alpha) * k_4) / (Ca(i) + k_4);
  217. }
  218.  
  219. double J_leak(double* f, int i)
  220. {
  221. return c_1 * v_2 * (c_0 / c_1 - (1 + 1 / c_1) * Ca(i));
  222. }
  223.  
  224. double J_pump(double* f, int i)
  225. {
  226. return v_3 * pow(Ca(i), 2) / (pow(k_3, 2) + pow(Ca(i), 2));
  227. }
  228.  
  229. double J_in(double* f, int i)
  230. {
  231. return v_5 + v_6 * pow(IP3(i), 2) / (pow(k_2, 2) + pow(IP3(i), 2));
  232. }
  233.  
  234. double J_out(double* f, int i)
  235. {
  236. return k_1 * Ca(i);
  237. }
  238.  
  239. double J_Glu(double* f, int i)
  240. {
  241. if (E_syn[i] == 0)
  242. {
  243. /*printf("J_Glu = %f\n", alpha_Glu / (1 + exp(-(G(i) - 0.4) / 0.01)));*/
  244. return alpha_Glu / (1 + exp(-(G(i) - 0.25) / 0.01));
  245. }
  246.  
  247. return 0;
  248. //return alpha_Glu / (1 + exp(-(G(i) - 0.25) / 0.01));
  249. }
  250.  
  251. double alpha_m(double* f, int i)
  252. {
  253. return 0.182 * (V(i) + 35) / (1 - exp(-(V(i) + 35) / 9));
  254. }
  255.  
  256. double beta_m(double* f, int i)
  257. {
  258. return -0.124 * (V(i) + 35) / (1 - exp((V(i) + 35) / 9));
  259. }
  260.  
  261. double alpha_n(double* f, int i)
  262. {
  263. return 0.02 * (V(i) - 25) / (1 - exp(-(V(i) - 25) / 9));
  264. }
  265.  
  266. double beta_n(double* f, int i)
  267. {
  268. return -0.002 * (V(i) - 25) / (1 - exp((V(i) - 25) / 9));
  269. }
  270.  
  271. double alpha_h(double* f, int i)
  272. {
  273. return 0.25 * exp(-(V(i) + 90) / 12);
  274. }
  275.  
  276. double beta_h(double* f, int i)
  277. {
  278. return 0.25 * exp((V(i) + 62) / 6) / exp((V(i) + 90) / 12);
  279. }
  280.  
  281. double UllahJung_HodgkinHuxley(int i, double* f, double t)
  282. {
  283. int in = i / Equations_per_node;
  284. int il = i % Equations_per_node;
  285.  
  286. switch (il)
  287. {
  288. case 0: // Ca
  289. {
  290. double sum_1 = 0;
  291.  
  292. /*for (int j = 0; j < Node_count; j++)
  293. {
  294. sum_1 += d_Ca * (Ca(j) - Ca(in));
  295. }*/
  296.  
  297. for (int j = 0; j < C_A[in]; j++)
  298. {
  299. sum_1 += d_Ca * (Ca((int)B_A[in][j]) - Ca(in));
  300. }
  301.  
  302. return J_channel(f, in) - J_pump(f, in) + J_leak(f, in) + J_in(f, in) - J_out(f, in) + sum_1;
  303. }
  304.  
  305. case 1: // IP3
  306. {
  307. double sum_2 = 0;
  308.  
  309. /*for (int j = 0; j < Node_count; j++)
  310. {
  311. sum_2 += d_IP3 * (IP3(j) - IP3(in));
  312. }*/
  313.  
  314. for (int j = 0; j < C_A[in]; j++)
  315. {
  316. sum_2 += d_IP3 * (IP3((int)B_A[in][j]) - IP3(in));
  317. }
  318.  
  319. return (IP3_star - IP3(in)) / tau_IP3 + J_PLC(f, in) + sum_2 + J_Glu(f, in);
  320. }
  321.  
  322. case 2: // z
  323. {
  324. return a_2 * (d_2 * (IP3(in) + d_1) / (IP3(in) + d_3) * (1 - z(in)) - Ca(in) * z(in));
  325. }
  326.  
  327. case 3: // G
  328. {
  329. return -alpha_G * G(in) + beta_G * (1 / (1 + exp(-V(in) / 0.5)));
  330. }
  331.  
  332. case 4: // V
  333. {
  334. double I_syn = 0;
  335.  
  336. /*for (int j = 0; j < Node_count; j++)
  337. {
  338. //sum += A[in][j] * g_syn * (V(in) - V_old(j, tau[in][j]));
  339. //sum += A[in][j] * g_syn * (V(j) - V(in));
  340. //sum += A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old(j, tau[in][j]) / k_syn));
  341. //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
  342. //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
  343. I_syn += A_N[in][j] * g_syn * (E_syn[in] - V(in)) / (1 + exp(-(V(j) / k_syn)));
  344. //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);
  345. }*/
  346.  
  347. /*for (int j = 0; j < C[in]; j++)
  348. {
  349. sum += sigma[in][(int)B[in][j]] * (V((int)B[in][j]) - V(in));
  350. }*/
  351.  
  352. for (int j = 0; j < C_N[in]; j++)
  353. {
  354. if (Ca(in) >= 0.3)
  355. {
  356. 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)
  357. }
  358. //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
  359. //I_syn += g_syn * (1 + g_astro * Ca(in)) * (E_syn[in] - V(in)) / (1 + exp(-(V(j) / k_syn))); // версия без с V (без V_old)
  360. //I_syn += g_syn * (E_syn[in] - V(in)) / (1 + exp(-(V(j) / k_syn))); // версия без с V (без V_old), упрощенная версия
  361. // sum_3 += g_syn * (1 + g_astro * Ca(in)) * (E_syn[i] - V(i)) / (1 + exp(-(V(j) / k_syn))); // образец из старой версии
  362. else
  363. {
  364. 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)
  365. }
  366. //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));*/
  367. //I_syn += g_syn * (E_syn[(int)B_N[in][j]] - V(in)) / (1 + exp(-(V(j) / k_syn))); // версия с V (без V_old)
  368. //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), упрощенная версия !!!
  369. //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), упрощенная версия !!!
  370. //sum += g_syn * (V((int)B[in][j]) - V(in)); // устаревшая часть, нужна для проверки разностной схемы
  371. //sum += A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old((int)B[in][j]) / k_syn)); // устаревшая часть
  372. //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));
  373. //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));
  374. //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
  375. //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
  376. // i up, j down
  377. /*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]]);*/
  378. /*printf("i = %d\t V_old = %f\t exp = %f\n", in, Vold, ee);*/
  379. }
  380. /*printf("i = %d\t sum = %f\n", in, sum);*/
  381.  
  382. if (enable_I_syn_out)
  383. fprintf(fp_I_syn, i == Equations_count - 1 ? "%f" : "%f\t", I_syn);
  384.  
  385. 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(35, t)*/ + I_syn) / C_m); // V
  386. }
  387.  
  388. case 5: // m
  389. {
  390. return 1000 * (alpha_m(f, in) * (1 - m(in)) - beta_m(f, in) * m(in)); // m
  391. }
  392.  
  393. case 6: // n
  394. {
  395. return 1000 * (alpha_n(f, in) * (1 - n(in)) - beta_n(f, in) * n(in)); // n
  396. }
  397.  
  398. case 7: // h
  399. {
  400. return 1000 * (alpha_h(f, in) * (1 - h(in)) - beta_h(f, in) * h(in)); // h
  401. }
  402. }
  403.  
  404. return 0;
  405. }
  406.  
  407. void RungeKutta(double t, double dt, double* f, double* f_next)
  408. {
  409. double k[Equations_count][4];
  410. double phi_k1[Equations_count];
  411. double phi_k2[Equations_count];
  412. double phi_k3[Equations_count];
  413.  
  414. // k1
  415. #pragma omp parallel for
  416. for (int i = 0; i < Equations_count; i++)
  417. {
  418. k[i][0] = UllahJung_HodgkinHuxley(i, f, t) * dt;
  419. phi_k1[i] = f[i] + k[i][0] / 2;
  420. k[i][1] = UllahJung_HodgkinHuxley(i, phi_k1, t) * dt;
  421. phi_k2[i] = f[i] + k[i][1] / 2;
  422. k[i][2] = UllahJung_HodgkinHuxley(i, phi_k2, t) * dt;
  423. phi_k3[i] = f[i] + k[i][2] / 2;
  424. k[i][3] = UllahJung_HodgkinHuxley(i, phi_k3, t) * dt;
  425. f_next[i] = f[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
  426. }
  427.  
  428. //for (int i = 0; i < Equations_count; i++)
  429. // phi_k1[i] = f[i] + k[i][0] / 2;
  430.  
  431. // k2
  432. //for (int i = 0; i < Equations_count; i++)
  433. // k[i][1] = UllahJung_HodgkinHuxley(i, phi_k1, t) * dt;
  434.  
  435.  
  436. //for (int i = 0; i < Equations_count; i++)
  437. // phi_k2[i] = f[i] + k[i][1] / 2;
  438.  
  439. // k3
  440. //for (int i = 0; i < Equations_count; i++)
  441. // k[i][2] = UllahJung_HodgkinHuxley(i, phi_k2, t) * dt;
  442.  
  443.  
  444. //for (int i = 0; i < Equations_count; i++)
  445. // phi_k3[i] = f[i] + k[i][2] / 2;
  446.  
  447. //enable_I_syn_out = true;
  448.  
  449. // k4
  450. //for (int i = 0; i < Equations_count; i++)
  451. // k[i][3] = UllahJung_HodgkinHuxley(i, phi_k3, t) * dt;
  452.  
  453. //enable_I_syn_out = false;
  454.  
  455. //for (int i = 0; i < Equations_count; i++)
  456. // f_next[i] = f[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
  457. }
  458.  
  459. void CopyArray(double* source, double* target, int N)
  460. {
  461. for (int i = 0; i < N; i++)
  462. target[i] = source[i];
  463. }
  464.  
  465. bool Approximately(double a, double b)
  466. {
  467. if (a < 0)
  468. a = -a;
  469.  
  470. if (b < 0)
  471. b = -b;
  472.  
  473. return a - b <= 0.000001;
  474. }
  475.  
  476. //bool CheckSameLine(int i, int j)
  477. //{
  478. // return i / Node_wire_width == j / Node_wire_width;
  479. //}
  480. //
  481. //bool IsWireNeighbors(int i, int j)
  482. //{
  483. // if (CheckSameLine(i, j) && (i == j - 1 || i == j + 1))
  484. // return true;
  485. //
  486. // if (i == j - Node_wire_width || i == j + Node_wire_width)
  487. // return true;
  488. //
  489. // return false;
  490. //}
  491.  
  492. // http://preshing.com/20111007/how-to-generate-random-timings-for-a-poisson-process/
  493. double nextTime(double rateParameter)
  494. {
  495. return -log(1.0 - (double)rand() / (RAND_MAX)) / rateParameter;
  496. }
  497.  
  498. void GenerateRandomMeander(int i, double min_start_time)
  499. {
  500. double offset = nextTime(Freq);
  501.  
  502. if (offset < 0)
  503. {
  504. int a = 0;
  505. }
  506.  
  507. Meander_start_from_zero[i] = min_start_time + offset;
  508. Meander_width[i] = Duration;
  509. Meander_height[i] = RandomD(Min_magintude, Max_magintude);
  510. }
  511.  
  512. void FillAMatrixZero()
  513. {
  514. for (int i = 0; i < Node_count; i++)
  515. {
  516. for (int j = 0; j < Node_count; j++)
  517. {
  518. A_A[i][j] = 0;
  519. A_N[i][j] = 0;
  520. }
  521. }
  522. }
  523.  
  524. void FillAstrociteMatrix()
  525. {
  526. for (int i = 0; i < Node_count; i++)
  527. {
  528. for (int j = 0; j < Node_count; j++)
  529. {
  530. if (i == j)
  531. {
  532. A_A[i][j] = 0;
  533. continue;
  534. }
  535.  
  536. if (i > j)
  537. {
  538. A_A[i][j] = A_A[j][i];
  539. continue;
  540. }
  541.  
  542. if (i == 0 && j == Node_count - 1)
  543. {
  544. A_A[i][j] = 1;
  545. continue;
  546. }
  547.  
  548. if (i == Node_count - 1 && j == 0)
  549. {
  550. A_A[i][j] = 1;
  551. continue;
  552. }
  553.  
  554. if (i == j - 1 || i == j + 1)
  555. {
  556. A_A[i][j] = 1;
  557. continue;
  558. }
  559. }
  560. }
  561. //A_A[0][1] = 0; // only for debug. diffusion Ca test
  562. //A_A[1][0] = 0; // only for debug. diffusion Ca test
  563. //A_A[1][3] = 0;
  564. //A_A[3][1] = 0;
  565. //A_A[0][2] = 0;
  566. //A_A[2][0] = 0;
  567. }
  568.  
  569. void FillBCMatrix_A()
  570. {
  571. for (int i = 0; i < Node_count; i++)
  572. {
  573. int bIndex = 0;
  574. C_A[i] = 0;
  575. for (int j = 0; j < Node_count; j++)
  576. {
  577. if (A_A[i][j] == 1)
  578. {
  579. B_A[i][bIndex] = j;
  580. bIndex++;
  581. C_A[i]++;
  582. }
  583. }
  584. }
  585. }
  586.  
  587. void FillBCMatrix_N()
  588. {
  589. for (int i = 0; i < Node_count; i++)
  590. {
  591. int bIndex = 0;
  592. C_N[i] = 0;
  593. for (int j = 0; j < Node_count; j++)
  594. {
  595. if (A_N[i][j] == 1)
  596. {
  597. B_N[i][bIndex] = j;
  598. bIndex++;
  599. C_N[i]++;
  600. }
  601. }
  602. }
  603. }
  604.  
  605. bool IsWireNeighbors(int i, int j, int deep)
  606. {
  607. int j_border_left = j - deep < 0 ? j + Node_count : j;
  608. int j_border_right = j + deep >= Node_count ? j - Node_count : j;
  609.  
  610. if (i == j_border_left - deep || i == j_border_right + deep)
  611. {
  612. return true;
  613. }
  614.  
  615. return false;
  616. }
  617.  
  618. void FillNeuronMatrix()
  619. {
  620. for (int i = 0; i < Node_count; i++)
  621. {
  622. for (int j = 0; j < Node_count; j++)
  623. {
  624. if (i == j)
  625. {
  626. A_N[i][j] = 0;
  627. continue;
  628. }
  629.  
  630. if (i > j)
  631. {
  632. A_N[i][j] = A_N[j][i];
  633. continue;
  634. }
  635.  
  636. for (int deep = 1; deep <= MaxDeep; deep++)
  637. {
  638. if (IsWireNeighbors(i, j, deep))
  639. A_N[i][j] = 1;
  640. }
  641. }
  642. }
  643. }
  644.  
  645. void RandomizeNeuronMatrix()
  646. {
  647. //srand(time(NULL));
  648.  
  649. for (int i = 0; i < Node_count; i++)
  650. {
  651. //if (i == Node_count / 2)
  652. // srand(time(NULL));
  653.  
  654. for (int link = 0; link < MaxDeep * 2; link++)
  655. {
  656. double x = RandomD(0, 1);
  657.  
  658. if (x > p_rewir)
  659. continue;
  660.  
  661. int rndJ;
  662.  
  663. do
  664. {
  665. rndJ = RandomI(0, Node_count);
  666. } while (i == rndJ || A_N[i][rndJ] == 1);
  667.  
  668. int rndJ_last;
  669.  
  670. do
  671. {
  672. rndJ_last = RandomI(i - MaxDeep - 1, i + MaxDeep + 1);
  673.  
  674. if (rndJ_last < 0)
  675. rndJ_last += Node_count;
  676. else if (rndJ_last >= Node_count)
  677. rndJ_last -= Node_count;
  678.  
  679. } while (i == rndJ_last || A_N[i][rndJ_last] == 0);
  680.  
  681. A_N[i][rndJ_last] = 0;
  682.  
  683. A_N[i][rndJ] = 1;
  684. }
  685. }
  686. }
  687.  
  688. void FillVOldFromCurrent()
  689. {
  690. for (int i = 0; i < Node_count; i++)
  691. for (int j = 0; j < Max_delay; j++)
  692. V_old_array[i][j] = V(i);
  693. }
  694.  
  695. void UpdateVOld()
  696. {
  697. for (int i = 0; i < Node_count; i++)
  698. {
  699. for (int j = 1; j < Max_delay; j++)
  700. V_old_array[i][j - 1] = V_old_array[i][j];
  701.  
  702. V_old_array[i][Max_delay - 1] = V(i);
  703. }
  704. }
  705.  
  706. //void FillFullTauMatrix()
  707. //{
  708. // for (int i = 0; i < Node_count; i++)
  709. // {
  710. // for (int j = 0; j < Node_count; j++)
  711. // {
  712. // if (i < Node_count || j < Node_count)
  713. // {
  714. // tau[i][j] = 0;
  715. // continue;
  716. // }
  717. //
  718. // int i_neuron = i - Node_count;
  719. // int j_neuron = j - Node_count;
  720. //
  721. // int i_wire_x = i_neuron / Node_wire_width;
  722. // int i_wire_y = i_neuron % Node_wire_width;
  723. //
  724. // int j_wire_x = j_neuron / Node_wire_width;
  725. // int j_wire_y = j_neuron % Node_wire_width;
  726. //
  727. // double distance_max = sqrt(2.) * (Node_wire_width - 1);
  728. // 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));
  729. //
  730. // tau[i][j] = (tau_min + distance / (distance_max) * (tau_max - tau_min)) * ms_to_step;
  731. // }
  732. // }
  733. //}
  734. //
  735. //void FillTauMatrix()
  736. //{
  737. // for (int i = 0; i < Node_count; i++)
  738. // {
  739. // for (int j = 0; j < Node_count; j++)
  740. // {
  741. // if (i == j || A_N[i][j] == 0)
  742. // {
  743. // tau[i][j] = 0;
  744. // continue;
  745. // }
  746. //
  747. // int i_neuron = i;
  748. // int j_neuron = j;
  749. //
  750. // int i_wire_x = i_neuron / Node_wire_width;
  751. // int i_wire_y = i_neuron % Node_wire_width;
  752. //
  753. // int j_wire_x = j_neuron / Node_wire_width;
  754. // int j_wire_y = j_neuron % Node_wire_width;
  755. //
  756. // double distance_max = sqrt(2.) * (Node_wire_width - 1);
  757. // 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));
  758. //
  759. // double t = (distance - 1) / (distance_max - 1);
  760. // tau[i][j] = (tau_min + t * (tau_max - tau_min)) * ms_to_step;
  761. // }
  762. // }
  763. //}
  764.  
  765. int main(int argc, char* argv[])
  766. {
  767. sscanf(argv[1], "%lf", &p_rewir);
  768. FILE* fp_p_rewir;
  769. fp_p_rewir = fopen("p_rewir.txt", "w");
  770. fprintf(fp_p_rewir, "%f\t", p_rewir);
  771. fclose(fp_p_rewir);
  772. printf("p_rewir = %f\n", p_rewir);
  773.  
  774. sscanf(argv[2], "%lf", &p_inhib);
  775. FILE* fp_p_inhib;
  776. fp_p_inhib = fopen("p_inhib.txt", "w");
  777. fprintf(fp_p_inhib, "%f\t", p_inhib);
  778. fclose(fp_p_inhib);
  779. printf("p_inhib = %f\n", p_inhib);
  780.  
  781. sscanf(argv[3], "%lf", &g_syn);
  782. FILE* fp_g_syn;
  783. fp_g_syn = fopen("g_syn.txt", "w");
  784. fprintf(fp_g_syn, "%f\t", g_syn);
  785. fclose(fp_g_syn);
  786.  
  787. double g_syn_real;
  788. g_syn_real = /*1 / (0.2 * Node_count) */ g_syn;
  789. printf("g_syn_real = %f\n", g_syn_real);
  790.  
  791. sscanf(argv[4], "%lf", &I_app_min);
  792. sscanf(argv[5], "%lf", &I_app_max);
  793.  
  794. FILE* fp0;
  795. FILE* fp_I_stim;
  796. FILE* fp_Ca;
  797. FILE* fp_IP3;
  798. //FILE *fp_z;
  799. FILE* fp_G;
  800. FILE* fp_V;
  801. //FILE *fp_m;
  802. //FILE *fp_n;
  803. //FILE *fp_h;
  804. //FILE *fp_V_spikes;
  805.  
  806. FILE* fp_Esyn;
  807. //srand(time(NULL));
  808.  
  809. //for (int i = 0; i < Node_count; i++)
  810. // V_old_length[i] = 0;
  811.  
  812. A_A = malloc(Node_count * sizeof(double));
  813. for (int i = 0; i < Node_count; i++)
  814. A_A[i] = malloc(Node_count * sizeof(double));
  815.  
  816. B_A = malloc(Node_count * sizeof(double));
  817. for (int i = 0; i < Node_count; i++)
  818. B_A[i] = malloc(Node_count * sizeof(double));
  819.  
  820. C_A = malloc(Node_count * sizeof(double));
  821.  
  822. A_N = malloc(Node_count * sizeof(double));
  823. for (int i = 0; i < Node_count; i++)
  824. A_N[i] = malloc(Node_count * sizeof(double));
  825.  
  826. B_N = malloc(Node_count * sizeof(double));
  827. for (int i = 0; i < Node_count; i++)
  828. B_N[i] = malloc(Node_count * sizeof(double));
  829.  
  830. C_N = malloc(Node_count * sizeof(double));
  831.  
  832. FillAMatrixZero();
  833. FillAstrociteMatrix();
  834. FillNeuronMatrix();
  835. RandomizeNeuronMatrix();
  836. FillBCMatrix_A();
  837. FillBCMatrix_N();
  838. //FillTauMatrix();
  839.  
  840. fp0 = fopen("A_A.txt", "w+");
  841. for (int i = 0; i < Node_count; i++)
  842. {
  843. for (int j = 0; j < Node_count; j++)
  844. {
  845. fprintf(fp0, "%d\t", (int)A_A[i][j]);
  846. }
  847. fprintf(fp0, "\n");
  848. }
  849. fclose(fp0);
  850.  
  851. fp0 = fopen("A_N.txt", "w+");
  852. for (int i = 0; i < Node_count; i++)
  853. {
  854. for (int j = 0; j < Node_count; j++)
  855. {
  856. fprintf(fp0, "%d\t", (int)A_N[i][j]);
  857. }
  858. fprintf(fp0, "\n");
  859. }
  860. fclose(fp0);
  861.  
  862. fp0 = fopen("tau.txt", "w+");
  863. for (int i = 0; i < Node_count; i++)
  864. {
  865. for (int j = 0; j < Node_count; j++)
  866. {
  867. fprintf(fp0, "%f\t", tau[i][j] / ms_to_step);
  868. }
  869. fprintf(fp0, "\n");
  870. }
  871. fclose(fp0);
  872.  
  873. //Пишем в файл число связей у каждого нейрона
  874. fp0 = fopen("links.txt", "w+");
  875. for (int i = 0; i < Node_count; i++)
  876. {
  877. int links_count = 0;
  878. for (int j = 0; j < Node_count; j++)
  879. {
  880. if (A_N[i][j] == 1)
  881. {
  882. links_count++;
  883. }
  884. }
  885. fprintf(fp0, "%d\n", (int)links_count);
  886. }
  887. fclose(fp0);
  888.  
  889. //setlocale(LC_NUMERIC, "French_Canada.1252");
  890. fp0 = fopen("test_Poisson.txt", "w+");
  891. for (int i = 0; i < 1000; i++)
  892. fprintf(fp0, "%f\n", nextTime(Freq));
  893. fclose(fp0);
  894.  
  895. fp0 = fopen("B_A.txt", "w+");
  896. for (int i = 0; i < Node_count; i++)
  897. {
  898. for (int j = 0; j < C_A[i]; j++)
  899. {
  900. fprintf(fp0, "%d\t", (int)B_A[i][j]);
  901. }
  902. fprintf(fp0, "\n");
  903. }
  904. fclose(fp0);
  905.  
  906. fp0 = fopen("B_N.txt", "w+");
  907. for (int i = 0; i < Node_count; i++)
  908. {
  909. for (int j = 0; j < C_N[i]; j++)
  910. {
  911. fprintf(fp0, "%d\t", (int)B_N[i][j]);
  912. }
  913. fprintf(fp0, "\n");
  914. }
  915. fclose(fp0);
  916.  
  917. fp0 = fopen("C_A.txt", "w+");
  918. for (int i = 0; i < Node_count; i++)
  919. {
  920. fprintf(fp0, "%d\n", (int)C_A[i]);
  921. }
  922. fclose(fp0);
  923.  
  924. fp0 = fopen("C_N.txt", "w+");
  925. for (int i = 0; i < Node_count; i++)
  926. {
  927. fprintf(fp0, "%d\n", (int)C_N[i]);
  928. }
  929. fclose(fp0);
  930.  
  931. /*for (int i = 0; i < 6; i++)
  932. {
  933. v_4[i] = 0.6;
  934. }*/
  935. for (int i = 0; i < Node_count; i++)
  936. {
  937. v_4[i] = 0.4; // 0.4
  938. }
  939.  
  940. // Initial values
  941. /*for (int i = 0; i < Equations_count; i++)
  942. {
  943. f[i] = 0;
  944. }*/
  945.  
  946. FILE* fp_I_app;
  947. fp_I_app = fopen("I_app.txt", "w");
  948. for (int i = 0; i < Node_count; i++)
  949. {
  950. I_app[i] = RandomD(I_app_min, I_app_max);
  951. fprintf(fp_I_app, "%f\n", I_app[i]);
  952. }
  953. fclose(fp_I_app);
  954.  
  955. for (int i = 0; i < Node_count; i++) // init array for all nodes
  956. {
  957. SetG(i, 0); // G
  958. }
  959.  
  960. double percent_stable_state = 0.4; // 0.40
  961. double eps_persent = 0.05; //0.05
  962.  
  963. double Ca0 = 0.07;
  964. double IP30 = 0.16;
  965. double z0 = 0.67;
  966.  
  967. // Initial values at t = 0
  968. /*for (int i = 0; i < Node_count; i++)
  969. {
  970. SetCa(i, Ca0); // Ca
  971. SetIP3(i, IP30); // IP3
  972. Setz(i, z0); // z
  973. }*/
  974. for (int i = 0; i < Node_count; i++)
  975. {
  976. /*SetCa(i, Ca0 + RandomD(-Ca0 * eps_persent, Ca0 * eps_persent)); // Ca
  977. SetIP3(i, IP30 + RandomD(-IP30 * eps_persent, IP30 * eps_persent)); // IP3
  978. Setz(i, z0 + RandomD(-z0 * eps_persent, z0 * eps_persent)); // z */
  979. SetCa(i, Ca0); // Ca
  980. SetIP3(i, IP30); // IP3
  981. Setz(i, z0); // z
  982. }
  983.  
  984. /*for (int i = 0; i < Node_count; i++) // init array for all nodes
  985. {
  986. SetV(i, V1); // V
  987. Setm(i, m1); // m
  988. Setn(i, n1); // n
  989. Seth(i, h1); // h
  990. }*/
  991.  
  992. double V0 = -58.7085;
  993. double m0 = 0.0953;
  994. double n0 = 0.000913;
  995. double h0 = 0.3662;
  996.  
  997. double V1 = 14.8409;
  998. double m1 = 0.9174;
  999. double n1 = 0.0140;
  1000. double h1 = 0.0539;
  1001.  
  1002. /*for (int i = 0; i < Node_count; i++) // init only for neurons
  1003. {
  1004. double random = RandomD(0, 1);
  1005.  
  1006. SetV(i, random < percent_stable_state ? V0 + RandomD(-V0 * eps_persent, V0 * eps_persent) : V1 + RandomD(-V1 * eps_persent, V1 * eps_persent)); // V
  1007. Setm(i, random < percent_stable_state ? m0 + RandomD(-m0 * eps_persent, m0 * eps_persent) : m1 + RandomD(-m1 * eps_persent, m1 * eps_persent)); // m
  1008. Setn(i, random < percent_stable_state ? n0 + RandomD(-n0 * eps_persent, n0 * eps_persent) : n1 + RandomD(-n1 * eps_persent, n1 * eps_persent)); // n
  1009. Seth(i, random < percent_stable_state ? h0 + RandomD(-h0 * eps_persent, h0 * eps_persent) : h1 + RandomD(-h1 * eps_persent, h1 * eps_persent)); // h
  1010. }*/
  1011.  
  1012. for (int i = 0; i < Node_count; i++) // init only for neurons
  1013. {
  1014. /*SetV(i, RandomD(-80, 20)); // V
  1015. Setm(i, RandomD(0, 1)); // m
  1016. Setn(i, RandomD(0, 1)); // n
  1017. Seth(i, RandomD(0, 1)); // h*/
  1018. SetV(i, V0); // V
  1019. Setm(i, m0); // m
  1020. Setn(i, n0); // n
  1021. Seth(i, h0); // h
  1022. }
  1023.  
  1024. double E_syn0 = 0; // Excitatory neuron
  1025. double E_syn1 = -90; // Inhibitory neuron
  1026.  
  1027. fp_Esyn = fopen("results_E_syn.txt", "w+");
  1028. for (int i = 0; i < Node_count; i++)
  1029. {
  1030. E_syn[i] = E_syn0;
  1031.  
  1032. double x = RandomD(0, 1);
  1033.  
  1034. if (x > p_inhib)
  1035. {
  1036. fprintf(fp_Esyn, "%f\n", E_syn[i]);
  1037. continue;
  1038. }
  1039.  
  1040. E_syn[i] = E_syn1;
  1041.  
  1042. fprintf(fp_Esyn, "%f\n", E_syn[i]);
  1043. }
  1044. fclose(fp_Esyn);
  1045.  
  1046. for (int i = 0; i < Node_count; i++)
  1047. {
  1048. GenerateRandomMeander(i, 0);
  1049. last_meander_end[i] = Meander_start_from_zero[i] + Duration;
  1050. }
  1051.  
  1052. const double t_start = 0;
  1053. const double t_max = 40.0; // 100 msec = 0.1 sec // 240 // 270
  1054. const double dt = 0.00005; // 0.01 msec = 0.00001 sec; 0.1 msec = 0.0001 sec; 1 msec = 0.001 sec // 0.000025
  1055.  
  1056. double t = t_start;
  1057.  
  1058. //fp0 = fopen("results.txt", "w+");
  1059. //setlocale(LC_NUMERIC, "French_Canada.1252");
  1060.  
  1061. clock_t start_rk4, end_rk4;
  1062. start_rk4 = clock();
  1063. int lastPercent = -1;
  1064.  
  1065. FillVOldFromCurrent();
  1066.  
  1067. fp_I_stim = fopen("results_I_stim.txt", "w+");
  1068. fp_I_syn = fopen("results_I_syn.txt", "w+");
  1069. fp_Ca = fopen("results_Ca.txt", "w+");
  1070. fp_IP3 = fopen("results_IP3.txt", "w+");
  1071. //fp_z = fopen("results_z.txt", "w+");
  1072. fp_G = fopen("results_G.txt", "w+");
  1073. fp_V = fopen("results_V.txt", "w+");
  1074. //fp_m = fopen("results_m.txt", "w+");
  1075. //fp_n = fopen("results_n.txt", "w+");
  1076. //fp_h = fopen("results_h.txt", "w+");
  1077. //fp_V_spikes = fopen("results_V_spikes.txt", "w+");
  1078.  
  1079. while (t < t_max || Approximately(t, t_max))
  1080. {
  1081. fprintf(fp_I_stim, "%f\t", t);
  1082. fprintf(fp_Ca, "%f\t", t);
  1083. fprintf(fp_IP3, "%f\t", t);
  1084. //fprintf(fp_z, "%f\t", t);
  1085. fprintf(fp_G, "%f\t", t);
  1086. fprintf(fp_V, "%f\t", t);
  1087. //fprintf(fp_m, "%f\t", t);
  1088. //fprintf(fp_n, "%f\t", t);
  1089. //fprintf(fp_h, "%f\t", t);
  1090. //fprintf(fp_V_spikes, "%f\t", t);
  1091.  
  1092. for (int i = 0; i < Node_count; i++)
  1093. {
  1094. if (t > last_meander_end[i])
  1095. {
  1096. GenerateRandomMeander(i, t);
  1097. last_meander_end[i] = Meander_start_from_zero[i] + Duration;
  1098. }
  1099.  
  1100. fprintf(fp_I_stim, "%f\t", I_stim(i, t));
  1101. }
  1102. fprintf(fp_I_stim, "\n");
  1103.  
  1104. for (int i = 0; i < Equations_count; i += Equations_per_node)
  1105. fprintf(fp_Ca, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // Ca
  1106.  
  1107. for (int i = 1; i < Equations_count; i += Equations_per_node)
  1108. fprintf(fp_IP3, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // IP3
  1109.  
  1110. //for (int i = 2; i < Equations_count; i += Equations_per_node)
  1111. // fprintf(fp_z, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // z
  1112.  
  1113. for (int i = 3; i < Equations_count; i += Equations_per_node)
  1114. fprintf(fp_G, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // G
  1115.  
  1116. for (int i = 4; i < Equations_count; i += Equations_per_node)
  1117. fprintf(fp_V, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // V
  1118.  
  1119. //for (int i = 5; i < Equations_count; i += Equations_per_node)
  1120. // fprintf(fp_m, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // m
  1121.  
  1122. //for (int i = 6; i < Equations_count; i += Equations_per_node)
  1123. // fprintf(fp_n, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // n
  1124.  
  1125. //for (int i = 7; i < Equations_count; i += Equations_per_node)
  1126. // fprintf(fp_h, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // h
  1127.  
  1128. fprintf(fp_Ca, "\n");
  1129. fprintf(fp_IP3, "\n");
  1130. //fprintf(fp_z, "\n");
  1131. fprintf(fp_G, "\n");
  1132. fprintf(fp_V, "\n");
  1133. //fprintf(fp_m, "\n");
  1134. //fprintf(fp_n, "\n");
  1135. //fprintf(fp_h, "\n");
  1136.  
  1137.  
  1138. double f_next[Equations_count];
  1139.  
  1140. RungeKutta(t, dt, f, f_next);
  1141.  
  1142. /*for (int i = 0; i < Equations_count; i += Equations_per_node)
  1143. {
  1144. double diff = f_next[i] - f[i];
  1145.  
  1146. fprintf(fp_V_spikes, i == Equations_count - 1 ? "%d" : "%d\t", diff < 0 && f_diff[i] > 0 && f[i] > -10 ? 1 : 0);
  1147.  
  1148. f_diff[i] = diff;
  1149. }*/
  1150.  
  1151. //fprintf(fp_V_spikes, "\n");
  1152.  
  1153. CopyArray(f_next, f, Equations_count);
  1154.  
  1155. t += dt;
  1156.  
  1157. int percent = (int)(100 * (t - t_start) / (t_max - t_start));
  1158. if (percent != lastPercent)
  1159. {
  1160. printf("Progress: %d%%\n", percent);
  1161. lastPercent = percent;
  1162. }
  1163.  
  1164. //printf("V(24) = %f\t V_old(24) = %f\n", f[24*4], V_old(24));
  1165. UpdateVOld();
  1166.  
  1167. fprintf(fp_I_syn, "\n");
  1168. }
  1169.  
  1170. fclose(fp_I_stim);
  1171. fclose(fp_I_syn);
  1172. fclose(fp_Ca);
  1173. fclose(fp_IP3);
  1174. //fclose(fp_z);
  1175. fclose(fp_G);
  1176. fclose(fp_V);
  1177. //fclose(fp_m);
  1178. //fclose(fp_n);
  1179. //fclose(fp_h);
  1180. //fclose(fp_V_spikes);
  1181.  
  1182. end_rk4 = clock();
  1183. double extime_rk4 = (double)(end_rk4 - start_rk4) / CLOCKS_PER_SEC;
  1184. int minutes = (int)extime_rk4 / 60;
  1185. int seconds = (int)extime_rk4 % 60;
  1186. printf("\nExecution time is: %d minutes %d seconds\n ", minutes, seconds);
  1187.  
  1188. fp0 = fopen("time_exec.txt", "w+");
  1189. fprintf(fp0, "%f\n", extime_rk4);
  1190. fclose(fp0);
  1191.  
  1192. for (int i = 0; i < Node_count; i++)
  1193. free(A_A[i]);
  1194. free(A_A);
  1195.  
  1196. for (int i = 0; i < Node_count; i++)
  1197. free(B_A[i]);
  1198. free(B_A);
  1199.  
  1200. free(C_A);
  1201.  
  1202. for (int i = 0; i < Node_count; i++)
  1203. free(A_N[i]);
  1204. free(A_N);
  1205.  
  1206. for (int i = 0; i < Node_count; i++)
  1207. free(B_N[i]);
  1208. free(B_N);
  1209.  
  1210. free(C_N);
  1211. }
Advertisement
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!