API tools faq
paste
SpaceQuester

Untitled

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