Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #define _USE_MATH_DEFINES
- #include "math.h"
- #include <stdlib.h>
- #include <stdio.h>
- #include <locale.h>
- #include <time.h>
- #include <stdbool.h>
- #define N 4
- #define V f[0]
- #define m f[1]
- #define n f[2]
- #define h f[3]
- double f[N];
- const double C = 1;
- double g_K = 36;
- double g_Na = 120;
- double g_L = 0.3;
- double E_K = -77; // -77 // -12
- double E_Na = 55; // 55 // 115
- double E_L = -54.4; // -54.4 // 10
- double I_app = 8.65;
- bool Random_meander = true;
- const double Min_freq = 50;
- const double Max_freq = 500;
- const double Min_magintude = -1.5;
- const double Max_magintude = 1.5;
- const double Duration = 0.001;
- double Meander_start_from_zero = 5;
- double Meander_width = 1;
- double Meander_height = 2;
- double Meander_interval = 3;
- int RandomI(int min, int max)
- {
- return ((double)rand() / (RAND_MAX - 1)) * (max - min) + min;
- }
- double RandomD(double min, double max)
- {
- return ((double)rand() / RAND_MAX) * (max - min) + min;
- }
- double Min(double a, double b)
- {
- return a > b ? b : a;
- }
- double Max(double a, double b)
- {
- return a > b ? a : b;
- }
- double alpha_n(double f[N])
- {
- return 0.01 * (V + 55) / (1 - exp(-(V + 55) / 10));
- //return (10 - V) / (100 * (exp((10 - V) / 10)) - 1);
- }
- double beta_n(double f[N])
- {
- return 0.125 * exp(-(V + 65) / 80);
- //return 0.125 * exp(-V / 80);
- }
- double alpha_m(double f[N])
- {
- return 0.1 * (V + 40) / (1 - exp(-(V + 40) / 10));
- //return (25 - V) / (10 * (exp((25 - V) / 10) - 1));
- }
- double beta_m(double f[N])
- {
- return 4 * exp(-(V + 65) / 18);
- //return 4 * exp(-V / 18);
- }
- double alpha_h(double f[N])
- {
- return 0.07 * exp(-(V + 65) / 20);
- //return 0.07 * exp(-V / 20);
- }
- double beta_h(double f[N])
- {
- return 1 / (exp(-(V + 35) / 10) + 1);
- //return 1 / (exp((30 - V) / 10) + 1);
- }
- double I_stim(double t)
- {
- if (t < Meander_start_from_zero)
- return 0;
- t -= Meander_start_from_zero;
- t = fmod(t, Meander_width + Meander_interval);
- return t < Meander_width ? Meander_height : 0;
- }
- double HodgkinHuxley(int i, double f[N], double t)
- {
- switch (i)
- {
- case 0:
- return (g_Na * m * m * m * h * (E_Na - V) + g_K * n * n * n * n * (E_K - V) + g_L * (E_L - V) + I_app + I_stim(t)) / C;
- case 1:
- return alpha_m(f) * (1 - m) - beta_m(f) * m;
- case 2:
- return alpha_n(f) * (1 - n) - beta_n(f) * n;
- case 3:
- return alpha_h(f) * (1 - h) - beta_h(f) * h;
- }
- return 0;
- }
- void RungeKutta(double t, double dt, double f[N], double f_next[N])
- {
- double k[N][4];
- // k1
- for (int i = 0; i < N; i++)
- k[i][0] = HodgkinHuxley(i, f, t) * dt;
- double phi_k1[N];
- for (int i = 0; i < N; i++)
- phi_k1[i] = f[i] + k[i][0] / 2;
- // k2
- for (int i = 0; i < N; i++)
- k[i][1] = HodgkinHuxley(i, phi_k1, t) * dt;
- double phi_k2[N];
- for (int i = 0; i < N; i++)
- phi_k2[i] = f[i] + k[i][1] / 2;
- // k3
- for (int i = 0; i < N; i++)
- k[i][2] = HodgkinHuxley(i, phi_k2, t) * dt;
- double phi_k3[N];
- for (int i = 0; i < N; i++)
- phi_k3[i] = f[i] + k[i][2] / 2;
- // k4
- for (int i = 0; i < N; i++)
- k[i][3] = HodgkinHuxley(i, phi_k3, t) * dt;
- for (int i = 0; i < N; i++)
- f_next[i] = f[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
- }
- void CopyArray(double source[N], double target[N])
- {
- for (int i = 0; i < N; i++)
- target[i] = source[i];
- }
- bool Approximately(double a, double b)
- {
- if (a < 0)
- a = -a;
- if (b < 0)
- b = -b;
- return a - b <= 0.000001;
- }
- void GenerateRandomMeander(double min_start_time)
- {
- Meander_start_from_zero = min_start_time + RandomD(1. / Max_freq, 1. / Min_freq);
- Meander_width = Duration;
- Meander_height = RandomD(Min_magintude, Max_magintude);
- }
- int main(int argc, char *argv[])
- {
- FILE *fp0;
- srand(time(NULL));
- for (int i = 0; i < N; i++)
- f[i] = 0;
- const double t_start = 0;
- const double t_max = 100; //2000
- const double dt = 0.0001; //0.05
- double t = t_start;
- fp0 = fopen("results.txt", "w+");
- //setlocale(LC_NUMERIC, "French_Canada.1252");
- clock_t start_rk4, end_rk4;
- start_rk4 = clock();
- int lastPercent = -1;
- double last_meander_end;
- if (Random_meander)
- {
- GenerateRandomMeander(0);
- last_meander_end = Meander_start_from_zero + Duration;
- }
- while (t < t_max || Approximately(t, t_max))
- {
- if (Random_meander && t > last_meander_end)
- {
- GenerateRandomMeander(t);
- last_meander_end = Meander_start_from_zero + Duration;
- }
- fprintf(fp0, "%f\t", t);
- fprintf(fp0, "%f\t", I_stim(t));
- for (int i = 0; i < N; i++)
- fprintf(fp0, i == N - 1 ? "%f" : "%f\t", f[i]);
- fprintf(fp0, "\n");
- double f_next[N];
- RungeKutta(t, dt, f, f_next);
- CopyArray(f_next, f);
- t += dt;
- int percent = (int)(100 * (t - t_start) / (t_max - t_start));
- if (percent != lastPercent)
- {
- printf("Progress: %d%%\n", percent);
- lastPercent = percent;
- }
- }
- end_rk4 = clock();
- double extime_rk4 = (double)(end_rk4 - start_rk4) / CLOCKS_PER_SEC;
- int minutes = (int)extime_rk4 / 60;
- int seconds = (int)extime_rk4 % 60;
- printf("\nExecution time is: %d minutes %d seconds\n ", minutes, seconds);
- fclose(fp0);
- fp0 = fopen("time_exec.txt", "w+");
- fprintf(fp0, "%f\n", extime_rk4);
- fclose(fp0);
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement