https://github.com/LKANG777/Beta-Oscillation
Revision dee52d24f28de078d94272131840b51c5f336513 authored by Ling KANG on 10 March 2023, 03:30:58 UTC, committed by GitHub on 10 March 2023, 03:30:58 UTC
1 parent c0feac8
Tip revision: dee52d24f28de078d94272131840b51c5f336513 authored by Ling KANG on 10 March 2023, 03:30:58 UTC
Create LICENSE
Create LICENSE
Tip revision: dee52d2
simulation.cpp
//
// simulation.cpp
// ==============
// Simulation of 2D model of motor cortex
// consisting of local E-I modules that are
// coupled with distance-dependent connection
// strengths, and subject to fluctuating external
// inputs. Individual modules are described by
// an adaptive rate model. For details, please see
// the corresponding publication:
// Ling, Ranft & Hakim, eLife (2023)
//
//
// Source code created by KANGLING on 2021/11/15.
//
// Imports
#include<iostream>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include<random>
#include<algorithm>
// Initialize random number generator
std::random_device rd; int seed = rd();
std::mt19937 rng(seed);
std::uniform_real_distribution<double> real_random(-0.1,0.1);
std::normal_distribution<double> normal_distribution(0,1.0);
/* Structure definitions */
/*-----------------------*/
/*
Two-dimensional network structure
t_tau: distance-dependent propagation delay
weight: long-range excitatory connection strength
color: different kinds of modules
*/
typedef struct nw
{
int **t_tau;
double **weight;
int *color;
} snw;
/*
Rate model f-I curve and adaptive time scale
N_data: length of the data
fi_curve: [0] current; [1] firing rate; [2]: time scale
*/
typedef struct f_tau
{
int N_data;
double **fi_curve;
} sf_tau;
/* Function declarations */
/*-----------------------*/
/* Get current, firing rate and time scale from the tabulated data for adaptive
* rate model */
sf_tau *read_data(double i_min, double i_max,double d_i);
double function_get_frequency(int N,double current,double**fi_curve);
double function_get_tau(int N,double current,double**fi_curve);
double function_get_current(int N,double r,double**fi_curve);
/* Instantiate the network */
snw*network_2D(int N,double lambda,int tau,int sur_width );
/* Compute the external current inputs */
double function_i_e_ext(double i_e_ext0, double sigma_ext_e, double eta_ind);
double function_i_i_ext(double i_i_ext0, double sigma_ext_i, double eta_ind);
double function_eta(double eta,double tau_ext,double xi_ind,double xi_all,double eta_c );
/* The kinetic kernels for the synaptic currents */
double function_current_arise(int i, int NN, double current_arise,double tau_arise, double r);
double function_current_decay(int i, int NN, double current_decay, double tau_decay,double current_arise);
/* E-I rate model related functions */
double function_i_e(int i,int NN,double i_e, double tau_e, double i_e_ext, double omega_ee,double omega_ei, double *current_decay_e, double current_decay_i, double**weight);
double function_i_i(int i,int NN,double i_i, double tau_i, double i_i_ext, double omega_ie,double omega_ii, double *current_decay_e, double current_decay_i, double **weight);
double delay_function_i_e(int i,int NN,double i_e, double tau_e, double i_e_ext, double omega_ee,double omega_ei, double *current_decay_e, double current_decay_i,double**weight,double**current_decay_e_tem,int **t_tau);
double delay_function_i_i(int i,int NN,double i_i, double tau_i, double i_i_ext, double omega_ie,double omega_ii, double *current_decay_e, double current_decay_i,double **weight,double **current_decay_e_tem,int **t_tau);
/* Simulation */
double calculate(int max_step,int N_data,double**fi_curve,int NN,double **weight,double lambda,double i_e_ext0, double i_i_ext0,double tau_ext,double sigma_ext_e,double sigma_ext_i,double omega_ee, double omega_ei,double omega_ie,double omega_ii,int N_noise,int **t_tau,int max_t_tau,double tau,double eta_c, int *color );
/* Simulation */
int main()
{
/*Recording simulation time*/
clock_t startTime,endTime;
startTime = clock();
/* Generating network
N: network size
sur_width: width of the fixed boundary
tau: basic propagation delay (unit: ms/step)
lambda: excitatory connectivity range
weight: long-range excitatory connection strength
t_tau: distance-dependent propagation delay
color: label to distinguish freely evolving modules from frozen modules at the fixed boundary
*/
int N = 28; // total grid length
int sur_width = 2; // width of surrounding layer (fixed rate)
double tau = 130; // delay in units of time step h (h = 10us)
double lambda = 2.0; //
int NN = N*N; // total number of modules
snw* snwk;
double **weight;
int **t_tau;
int *color;
printf("Network parameters:\nN = %d (network size)\nsur_width = %d (width of fixed-rate module boundary)\nD = %.2f ms (delay between neighboring modules)\ntau = %.0f (D/simulation step,simulation step h=0.01)\nlambda = %.2f (exc. connectivity range)\n", N, sur_width, tau/100, tau, lambda);
// Building the network
snwk = network_2D(N,lambda,tau,sur_width );
weight = (snwk->weight);
t_tau = (snwk->t_tau);
color = (snwk->color);
// Writing the coupling strength and delay step between each two nodes to file:
FILE*fp_weight;
char filename_weight[256];
sprintf(filename_weight,"%d_%.2f_weight.txt",N,lambda);
fp_weight=fopen(filename_weight,"w");
for(int i=0;i<NN;i++)
{
double add_weight=0.0;
for(int j=0;j<NN;j++)
{
add_weight+=weight[i][j];
fprintf(fp_weight,"%f,",weight[i][j]);
}
fprintf(fp_weight,"%f,\n",add_weight);
}
// Determine max. value of delay that will be needed to be taken into account:
int max_t_tau=0;
for(int i=0;i<NN;i++)
{
for(int j=0;j<NN;j++)
{
if(t_tau[i][j]>max_t_tau) max_t_tau=t_tau[i][j];
fprintf(fp_weight,"%d,",t_tau[i][j]);
}
fprintf(fp_weight,"\n");
}
fclose(fp_weight);
printf("Max. number of past time steps needed for delay: max_t_tau=%d\n",max_t_tau);
// Record the color of the different modules/nodes (fixed nodes and the simulated nodes):
FILE * fp_color;
char filename_color[256];
sprintf(filename_color, "%d_%f_color.txt", N,lambda);
fp_color = fopen(filename_color, "w");
for (int i = 0; i < N; i++)
{
for (int j = 0; j <N; j++)
{
int jj= i*N+j;
fprintf(fp_color, "%d\t", color[jj]);
}
fprintf(fp_color, "\n");
}
fclose(fp_color);
/* Reading original tabulated data for rate model, i.e., f-I curve and adaptive time scale */
double i_min = -20.0;
double i_max = 100.0;
double d_i = 0.1;
sf_tau*sf_tauk;
double**fi_curve;
int N_data;
sf_tauk = read_data(i_min,i_max,d_i);
fi_curve = (sf_tauk->fi_curve);
N_data = (sf_tauk->N_data);
/* Estimated number of neurons per module, setting the strength of the finite-size noise due to Poissonian spiking*/
int N_noise = 20000;
printf("N_noise = %d (number of neurons per module setting finite-size noise strength)\n",N_noise);
/*Recurrent synaptic coupling strength (mV.s)*/
double omega_ee = 0.96;
double omega_ie = 1.0;
double omega_ei = 2.08;
double omega_ii = 0.87;
/*
External input
eta_c: proportion of global external inputs
nu (Hz): external input amplitude fluctuations
tau_ext (ms): correlation time of external input fluctuations
*/
double nu = 3;
double sigma_ext_e = nu*omega_ee;
double sigma_ext_i = 2*nu*omega_ie;
double eta_c = 0.4;
double tau_ext = 25;
/*Steady state*/
double r_e_s = 5.0;
double r_i_s = 10.0;
double i_e_s;
double i_i_s;
i_e_s = function_get_current(N_data, r_e_s,fi_curve);
i_i_s = function_get_current(N_data, r_i_s,fi_curve);
double i_e_ext0 = i_e_s-omega_ee*r_e_s+omega_ei*r_i_s;
double i_i_ext0 = i_i_s-omega_ie*r_e_s+omega_ii*r_i_s;
printf("\nConnection parameters:\nw_{ee} = %.2f\nw_{ei} = %.2f\nw_{ie} = %.2f\nw_{ii} = %.2f\n\nInput parameters:\nnu = %.2f\neta_c = %.2f\ntau_{ext} = %.2f\n", omega_ee,omega_ei,omega_ie,omega_ii,nu,eta_c ,tau_ext );
/* Setting up the simulation */
int duration = 100; // simulation duration in ms
int max_step = duration*100; // step size hardcoded to be 0.01 ms
printf("\nSimulation duration:\nT = %d ms (%d steps)\n\n", duration, max_step);
// The simulation...
calculate(max_step, N_data, fi_curve, NN, weight,lambda,i_e_ext0, i_i_ext0, tau_ext,sigma_ext_e,sigma_ext_i, omega_ee, omega_ei, omega_ie, omega_ii, N_noise,t_tau,max_t_tau,tau,eta_c ,color);
endTime = clock();
std::cout << "\nTotal time needed: " <<(double)(endTime - startTime) / CLOCKS_PER_SEC << "s" << std::endl;
// printf("end");
// getchar();
// return 0;
}
/*----------------------*/
/* Function definitions */
/*----------------------*/
/* Instatiating and initialising the network */
snw*network_2D(int N,double lambda,int tau,int sur_width )
{
int NN=N*N;
snw*snwk;
snwk=(snw*)malloc(1*sizeof(snw));
int **degree;
int **adjacent;
double **weight;
int **t_tau;
t_tau=(int**)malloc(NN*sizeof(int*));
degree=(int**)malloc(NN*sizeof(int*));
adjacent=(int**)malloc(NN*sizeof(int*));
weight=(double**)malloc(NN*sizeof(double*));
for (int i = 0; i < NN; i++)
{
adjacent[i] = (int*)malloc(NN * sizeof(int));
degree[i] = (int*)malloc(NN * sizeof(int));
weight[i] = (double*)malloc(NN * sizeof(double));
t_tau[i] = (int*)malloc(NN*sizeof(int));
for(int j=0;j<NN;j++)
{
weight[i][j]=0.0;
adjacent[i][j] = 0;
degree[i][j] = 0;
t_tau[i][j]=0;
}
}
for (int i = 0; i < NN; i++)
{
adjacent[i][i]=1;
degree[i][0]++;
degree[i][degree[i][0]] = i;
for (int j = i+1; j < NN; j++)
{
adjacent[i][j] = adjacent[j][i] = 1;
degree[i][0]++;
degree[i][degree[i][0]] = j;
degree[j][0]++;
degree[j][degree[j][0]] = i;
}
}
int *row;
int *col;
row = (int*)malloc(NN * sizeof(int));
col = (int*)malloc(NN * sizeof(int));
for (int i = 0; i < NN; i++)
{
row[i] = 0;
col[i] = 0;
}
for (int i = 0; i < NN; i++)
{
row[i] = int(i / N);
col[i] = i % N;
}
// Calculate the distance between units
double **distance;
distance = (double**)malloc(NN * sizeof(double*));
for (int i = 0; i < NN; i++)
{
distance[i] = (double*)malloc(NN * sizeof(double));
for (int j = 0; j < NN; j++)
{
distance[i][j] = 0;
}
}
int row_d;
int col_d;
for (int i = 0; i < NN; i++)
{
for (int j = 0; j <NN; j++)
{
if (abs(row[j] - row[i]) > N/2)
{
if (row[j] - row[i] > 0)
{
row_d = N - row[j] + row[i];
}
else
{
row_d = N - row[i] + row[j];
}
}
else
{
row_d = row[j] - row[i];
}
if (abs(col[j] - col[i]) > N/2)
{
if (col[j] - col[i] > 0)
{
col_d = N - col[j] + col[i];
}
else
{
col_d = N - col[i] + col[j];
}
}
else
{
col_d = col[j] - col[i];
}
distance[i][j] = sqrt(row_d * row_d + col_d * col_d);
// printf("%d,%d,%f\n", jj, i, distance[i][j]);
}
}
// Determine the distance-dependent weights between units
double add_weight=0;
for(int i=0;i<NN;i++)
{
for(int j=0;j<NN;j++)
{
double d=distance[i][j];
weight[i][j]=weight[j][i]= exp(-((d/lambda)*(d/lambda)));
}
if(i==0)
{
for(int j=0;j<NN;j++)
{
add_weight+=weight[i][j];
}
}
}
for(int i=0;i<NN;i++)
{
for(int j=0;j<NN;j++)
{
weight[i][j]= weight[i][j]/add_weight;
}
}
/*Finding the effective length of coupling*/
int R=N;
for(int i=0;i<1;i++)
{
for(int j=0;j<N;j++)
{
if (weight[i][j]<0.0000001)
{
R=j;
break;
}
}
}
// Determine the effective delay between connected modules
for(int i=0;i<NN;i++)
{
for(int j=0;j<NN;j++)
{
if (distance[i][j]<R)
{
t_tau[i][j]=t_tau[j][i]=(int)(tau*distance[i][j]);
}
}
}
// Write distances, weights, and delays to file
FILE * fp_net;
char filename_net[256];
sprintf(filename_net, "%d_%.2f_distance.txt", N,lambda);
fp_net = fopen(filename_net, "w");
for (int i = 0; i < NN; i++)
{
for (int j = 0; j <NN; j++)
{
fprintf(fp_net, "%d,%d,%f,%f,%d\n", i, j, distance[i][j]*distance[i][j],weight[i][j],t_tau[i][j]);
}
}
fclose(fp_net);
/*Distinguishing different kinds modules*/
int *color;
color=(int*)malloc(NN*sizeof(int));
for (int i = 0; i < NN; i++)
{
color[i] =0;
}
for (int i = 0; i < sur_width; i++)
{
for(int j=i*N; j<((i+1)*N);j++)
{
color[j]=1;
}
for(int j=0; j<N;j++)
{
int jj=j*N+i;
color[jj]=1;
}
}
for (int i = N-sur_width; i < N; i++)
{
for(int j=i*N; j<((i+1)*N);j++)
{
color[j]=1;
}
for(int j=0; j<N;j++)
{
int jj=j*N+i;
color[jj]=1;
}
}
(snwk->color) = color;
(snwk->t_tau) = t_tau;
(snwk->weight)=weight;
return snwk;
}
/* Read f-I curve and current-dependent timescale of rate model
* from tabulated data */
sf_tau *read_data(double i_min,double i_max,double di)
{
int N=(int)((i_max-i_min)/di)+1;
int n_bin=5; //reorder the bin;
int re_N=N*n_bin-n_bin+1;
sf_tau*sf_tauk;
sf_tauk = (sf_tau*)malloc(sizeof(sf_tau));
double**fi_curve;
double**upsampled_fi_curve;
fi_curve = (double**)malloc(N * sizeof(double*));
upsampled_fi_curve = (double**)malloc(re_N * sizeof(double*));
for (int i = 0; i < N; i++)
{
fi_curve[i] = (double*)malloc(3 * sizeof(double));
for (int j = 0; j < 3; j++)
{
fi_curve[i][j] = 0.0;
}
}
for (int i = 0; i < re_N; i++)
{
upsampled_fi_curve[i] = (double*)malloc(3 * sizeof(double));
for (int j = 0; j < 3; j++)
{
upsampled_fi_curve[i][j] = 0.0;
}
}
// Reading the tabulated data from file
FILE *fp;
fp = fopen("tabulated_fi_curve.txt", "r");
if (!fp)
{
printf("can't open file\n");
getchar();
exit(1);
}
for (int i = 0; i < N; i++)
{
for (int j = 0; j < 3; j++)
{
fscanf(fp, "%lf,", &fi_curve[i][j]);
// printf("%d,%f\n,", i,fi_curve[i][j]);
}
fi_curve[i][1]=fi_curve[i][1]*1000;
}
fclose(fp);
// Create an upsampled version
for (int i = 0; i < N; i++)
{
if(i<N-1)
{
for (int bin=0;bin<n_bin;bin++)
{
int re_i=i*n_bin+bin;
if(bin==0)
{
for (int j = 0; j < 3; j++)
{
upsampled_fi_curve[re_i][j]=fi_curve[i][j];
}
}
else
{
for (int j = 0; j < 3; j++)
{
upsampled_fi_curve[re_i][j]=(fi_curve[i+1][j]-fi_curve[i][j])*bin/n_bin+fi_curve[i][j];
}
}
}
}
}
upsampled_fi_curve[re_N-1][0]=fi_curve[N-1][0];
upsampled_fi_curve[re_N-1][1]=fi_curve[N-1][1];
upsampled_fi_curve[re_N-1][2]=fi_curve[N-1][2];
(sf_tauk->fi_curve) = upsampled_fi_curve;
(sf_tauk->N_data) = re_N;
return sf_tauk;
}
// Interpolate rate model timescale from data
double function_get_tau(int N,double current,double**fi_curve)
{
double dis_current=fi_curve[1][0]-fi_curve[0][0];
double div_bin=1000;
double unit_div_bin=dis_current/div_bin;
int dis=0;
double dis_tau=0.0;
double tau=0.0;
if(current<-20.0)
{
tau=fi_curve[0][2];
}
else if(current>100)
{
tau=fi_curve[N-1][2];
}
else{
for(int i=0;i<N;i++)
{
if((fi_curve[i][0]-current)>0.0000)
{
dis=(int)((current-fi_curve[i-1][0])/unit_div_bin);
dis_tau=fi_curve[i][2]-fi_curve[i-1][2];
tau=fi_curve[i-1][2]+dis_tau/div_bin*dis;
break;
}
}
}
return tau;
}
// Interpolate rate model firing rate from data
double function_get_frequency(int N,double current,double**fi_curve)
{
double dis_current=fi_curve[1][0]-fi_curve[0][0];
double div_bin=1000;
double unit_div_bin=dis_current/div_bin;
int dis=0;
double dis_r=0.0;
double r=0;
if(current<-20.0)
{
r=fi_curve[0][1];
}
else if(current>100)
{
r=fi_curve[N-1][1];
}
else{
for(int i=0;i<N;i++)
{
if((fi_curve[i][0]-current)>0.0000)
{
dis=(int)((current-fi_curve[i-1][0])/unit_div_bin);
dis_r=fi_curve[i][1]-fi_curve[i-1][1];
r=fi_curve[i-1][1]+dis_r/div_bin*dis;
break;
}
}
}
return r;
}
// Interpolate inverse f-I curve to get current related to given firing rate
double function_get_current(int N,double r,double**fi_curve)
{
double dis_r=fi_curve[1][1]-fi_curve[0][1];
double div_bin=1000;
double unit_div_bin=dis_r/div_bin;
int dis=0;
double dis_current=0.0;
double current=0;
if(r< 0)
{
r=fi_curve[0][0];
}
else if(r>262.2737)
{
r=fi_curve[N-1][0];
}
else{
for(int i=0;i<N;i++)
{
if((fi_curve[i][1]-r)>0.0000)
{
dis=(int)((r-fi_curve[i-1][1])/unit_div_bin);
dis_current=fi_curve[i][0]-fi_curve[i-1][0];
current=fi_curve[i-1][0]+dis_r/div_bin*dis;
break;
}
}
}
return current;
}
// External input to exc. populations
double function_i_e_ext(double i_e_ext0, double sigma_ext_e, double eta_ind )
{
double f_i_e_ext=0;
f_i_e_ext= i_e_ext0+sigma_ext_e*eta_ind;
return f_i_e_ext;
}
// External input to inh. populations
double function_i_i_ext(double i_i_ext0, double sigma_ext_i,double eta_ind)
{
double f_i_i_ext=0;
f_i_i_ext= i_i_ext0+sigma_ext_i*eta_ind;
return f_i_i_ext;
}
// Time derivative Ornstein-Uhlenbeck process (global + local contributions)
double function_eta(double eta,double tau_ext,double xi_ind,double xi_all,double eta_c )
{
double f_eta=0;
f_eta=1.0/tau_ext*(-eta+sqrt(tau_ext)*(sqrt(1.0-eta_c)*xi_ind+sqrt(eta_c)*xi_all)) ;
return f_eta;
}
double function_xi_ind(double xi_ind_pre,double tau_ext,double xi_ind,double xi_all,double eta_c )
{
double f_xi_ind=0;
f_xi_ind=1.0/tau_ext*(-xi_ind_pre+sqrt(tau_ext)*(sqrt(1.0-eta_c)*xi_ind)) ;
return f_xi_ind;
}
double function_xi_all(double xi_all_pre,double tau_ext,double xi_ind,double xi_all,double eta_c )
{
double f_xi_all=0;
f_xi_all=1.0/tau_ext*(-xi_all_pre+sqrt(tau_ext)*(sqrt(eta_c)*xi_all)) ;
return f_xi_all;
}
// Time derivative synaptic current (rise)
double function_current_arise( double current_arise,double tau_arise, double r)
{
double f_current_arise;
f_current_arise=1.0/tau_arise*(-current_arise+r);
return f_current_arise;
}
// Time derivative synaptic current (full kinetics)
double function_current_decay( double current_decay, double tau_decay,double current_arise)
{
double f_current_decay;
f_current_decay=1.0/tau_decay*(-current_decay+current_arise );
return f_current_decay;
}
// Time derivative driving current of exc. population (instantaneous)
double function_i_e(int i,int NN,double i_e, double tau_e, double i_e_ext, double omega_ee,double omega_ei, double *current_decay_e, double current_decay_i,double**weight)
{
double f_i_e=0;
double sum=0.0;
for (int j=0;j<NN;j++)
{
sum+=weight[i][j]*current_decay_e[j];
}
f_i_e=1.0/tau_e*(-i_e+i_e_ext+omega_ee*sum-omega_ei*current_decay_i);
return f_i_e;
}
// Time derivative driving current of inh. population (instantaneous)
double function_i_i(int i,int NN,double i_i, double tau_i, double i_i_ext, double omega_ie,double omega_ii, double *current_decay_e, double current_decay_i, double **weight)
{
double f_i_i=0;
double sum=0.0;
for (int j=0;j<NN;j++)
{
sum+=weight[i][j]*current_decay_e[j];
}
f_i_i=1.0/tau_i*(-i_i+i_i_ext+omega_ie*sum-omega_ii*current_decay_i);
return f_i_i;
}
// Time derivative driving current of exc. population (including delay)
double delay_function_i_e(int i,int NN,double i_e, double tau_e, double i_e_ext, double omega_ee,double omega_ei, double *current_decay_e, double current_decay_i,double**weight,double**current_decay_e_tem,int **t_tau)
{
double f_i_e=0;
double sum=0.0;
for (int j=0;j<NN;j++)
{
sum+=weight[i][j]*current_decay_e_tem[j][t_tau[i][j]];
}
f_i_e=1.0/tau_e*(-i_e+i_e_ext+omega_ee*sum-omega_ei*current_decay_i);
return f_i_e;
}
// Time derivative driving current of exc. population (including delay)
double delay_function_i_i(int i,int NN,double i_i, double tau_i, double i_i_ext, double omega_ie,double omega_ii, double *current_decay_e, double current_decay_i,double **weight,double **current_decay_e_tem,int **t_tau)
{
double f_i_i=0;
double sum=0.0;
for (int j=0;j<NN;j++)
{
sum+=weight[i][j]*current_decay_e_tem[j][t_tau[i][j]];
}
f_i_i=1.0/tau_i*(-i_i+i_i_ext+omega_ie*sum-omega_ii*current_decay_i);
return f_i_i;
}
// Actual simulation
double calculate(int max_step,int N_data,double**fi_curve,int NN,double **weight,double lambda,double i_e_ext0, double i_i_ext0,double tau_ext,double sigma_ext_e,double sigma_ext_i,double omega_ee, double omega_ei,double omega_ie,double omega_ii,int N_noise,int **t_tau,int max_t_tau,double tau,double eta_c, int *color )
{
/*Recording simulation time*/
clock_t startTime,endTime;
startTime = clock();
/*Recording simulation data*/
FILE*fp;
char filename[256];
sprintf(filename,"%d_%.2f_%d_%.2f_%.2f_%.2f_%.2f_%.2f_%.2f_%.2f_%.2f_current_tau_fix_e.txt", int(sqrt(NN)),tau*0.01,N_noise,lambda,omega_ee,omega_ei,omega_ie,omega_ii,sigma_ext_e/omega_ee,eta_c ,tau_ext);
fp=fopen(filename,"w");
FILE*fp_i;
char filename_i[256];
sprintf(filename_i,"%d_%.2f_%d_%.2f_%.2f_%.2f_%.2f_%.2f_%.2f_%.2f_%.2f_current_tau_fix_i.txt", int(sqrt(NN)),tau*0.01,N_noise,lambda,omega_ee,omega_ei,omega_ie,omega_ii,sigma_ext_e/omega_ee,eta_c ,tau_ext);
fp_i=fopen(filename_i,"w");
FILE*fp_xi_g;
char filename_xi_g[256];
sprintf(filename_xi_g,"%d_%.2f_%d_%.2f_%.2f_%.2f_%.2f_%.2f_%.2f_%.2f_%.2f_current_tau_fix_xi_g.txt", int(sqrt(NN)),tau*0.01,N_noise,lambda,omega_ee,omega_ei,omega_ie,omega_ii,sigma_ext_e/omega_ee,eta_c ,tau_ext);
fp_xi_g=fopen(filename_xi_g,"w");
// Initialize network rate model variables (currents, rates, timescales)
double *i_e,*i_i,*r_e,*r_i,*tau_e,*tau_i;
i_e=(double*)malloc(NN*sizeof(double));
i_i=(double*)malloc(NN*sizeof(double));
r_e=(double*)malloc(NN*sizeof(double));
tau_e=(double*)malloc(NN*sizeof(double));
r_i=(double*)malloc(NN*sizeof(double));
tau_i=(double*)malloc(NN*sizeof(double));
for(int i=0;i<NN;i++)
{
i_e[i]=-6.7;
i_i[i]=0;
r_e[i]=0.0;
tau_e[i]=0.0;
r_i[i]=0.0;
tau_i[i]=0.0;
}
double *i_e_pre,*i_i_pre;
i_e_pre=(double*)malloc(NN*sizeof(double));
i_i_pre=(double*)malloc(NN*sizeof(double));
for(int i=0;i<NN;i++)
{
i_e_pre[i]=0.0;
i_i_pre[i]=0.0;
}
/*External current*/
double *i_e_ext,*i_i_ext,*eta_ind,*xi_ind_pre,*eta_ind_pre,*xi_ind;
double xi_all_pre, xi_all;
i_e_ext=(double*)malloc(NN*sizeof(double));
i_i_ext=(double*)malloc(NN*sizeof(double));
//xi_all_pre=(double*)malloc(NN*sizeof(double));
eta_ind=(double*)malloc(NN*sizeof(double));
xi_ind_pre=(double*)malloc(NN*sizeof(double));
eta_ind_pre=(double*)malloc(NN*sizeof(double));
//xi_all=(double*)malloc(NN*sizeof(double));
xi_ind=(double*)malloc(NN*sizeof(double));
xi_all_pre=0;
for(int i=0;i<NN;i++)
{
i_e_ext[i]=0;
i_i_ext[i]=0;
//xi_all_pre[i]=0;
xi_all_pre=0;
xi_ind_pre[i]=real_random(rng);
//eta_ind[i]=xi_ind_pre[i]+xi_all_pre[0];
//eta_ind_pre[i]= xi_ind_pre[i]+xi_all_pre[0];
eta_ind[i]=xi_ind_pre[i]+xi_all_pre;
eta_ind_pre[i]= xi_ind_pre[i]+xi_all_pre;
//xi_all[i]=0;
xi_all=0;
xi_ind[i]=0;
}
/*Synatic current*/
double *current_decay_e,*current_arise_e,*current_decay_i,*current_arise_i;
double *current_decay_e_pre,*current_arise_e_pre,*current_decay_i_pre,*current_arise_i_pre;
current_decay_e=(double*)malloc(NN*sizeof(double));
current_arise_e=(double*)malloc(NN*sizeof(double));
current_decay_i=(double*)malloc(NN*sizeof(double));
current_arise_i=(double*)malloc(NN*sizeof(double));
current_decay_e_pre=(double*)malloc(NN*sizeof(double));
current_arise_e_pre=(double*)malloc(NN*sizeof(double));
current_decay_i_pre=(double*)malloc(NN*sizeof(double));
current_arise_i_pre=(double*)malloc(NN*sizeof(double));
for(int i=0;i<NN;i++)
{
current_decay_e[i]=real_random(rng);
current_arise_e[i]=real_random(rng);
current_decay_i[i]=real_random(rng);
current_arise_i[i]=real_random(rng);
current_decay_e_pre[i]=0;
current_arise_e_pre[i]=0;
current_decay_i_pre[i]=0;
current_arise_i_pre[i]=0;
}
double tau_decay=3.5,tau_arise=0.7,tau_lat=0.5;
double h=0.01;
double sqrt_h=sqrt(h);
/* Finite-size-noise of populations depends on number of neurons */
double N_e=0.8*N_noise;
double N_i=0.2*N_noise;
//** latency delay
int max_t_lat=int(tau_lat/h);
double **r_e_lat,**r_i_lat;
r_e_lat=(double**)malloc(NN*sizeof(double*));
r_i_lat=(double**)malloc(NN*sizeof(double*));
for (int i=0;i<NN;i++)
{
r_e_lat[i]=(double*)malloc((max_t_lat+1)*sizeof(double));
r_i_lat[i]=(double*)malloc((max_t_lat+1)*sizeof(double));
for (int j=0;j<(max_t_lat+1);j++)
{
r_e_lat[i][j]=0.0;
r_i_lat[i][j]=0.0;
}
}
int delay_step=max_t_tau+1;
/*Delay*/
double **current_decay_e_tem;
current_decay_e_tem=(double**)malloc(NN*sizeof(double*));
for (int i=0;i<NN;i++)
{
current_decay_e_tem[i]=(double*)malloc((max_t_tau+1)*sizeof(double));
for (int j=0;j<(max_t_tau+1);j++)
{
current_decay_e_tem[i][j]=0.0;
}
}
/*Simulation*/
for(int step=0;step<max_step+1;step++)
{
//**update current,frequence,tau,convey delay;
for(int i=0;i<NN;i++)
{
i_e_pre[i]=i_e[i];
i_i_pre[i]=i_i[i];
eta_ind_pre[i]=eta_ind[i];
current_arise_e_pre[i]=current_arise_e[i];
current_decay_e_pre[i]=current_decay_e[i];
current_arise_i_pre[i]=current_arise_i[i];
current_decay_i_pre[i]=current_decay_i[i];
// determine current-dependent deterministic firing rate
r_e[i]= function_get_frequency(N_data, i_e_pre[i],fi_curve);
// when finite-size noise is considered, sample spikes from
// Poisson distribution and determine stochastic firing rate
if(N_noise>0)
{
std::poisson_distribution<int> distribution_r_e(N_e*r_e[i]*h/1000);
r_e[i]=distribution_r_e(rng)/(N_e*h/1000);
}
tau_e[i]= function_get_tau(N_data, i_e_pre[i],fi_curve);
r_i[i]= function_get_frequency(N_data, i_i_pre[i],fi_curve);
if(N_noise>0)
{
std::poisson_distribution<int> distribution_r_i(N_i*r_i[i]*h/1000);
r_i[i]=distribution_r_i(rng)/(N_i*h/1000);
}
tau_i[i]= function_get_tau(N_data, i_i_pre[i],fi_curve);
r_e_lat[i][0]=r_e[i];
r_i_lat[i][0]=r_i[i];
}
//**calculate external current;
xi_all = normal_distribution(rng);
xi_all_pre = xi_all_pre + h*function_xi_all(xi_all_pre,tau_ext,xi_ind[0]/sqrt_h,xi_all/sqrt_h, eta_c );
for(int i=0;i<NN;i++)
{
//xi_ali]= normal_distribution(rng);
xi_ind[i]= normal_distribution(rng);
eta_ind[i]=eta_ind_pre[i]+h*function_eta(eta_ind_pre[i],tau_ext,xi_ind[i]/sqrt_h,xi_all/sqrt_h, eta_c );
xi_ind_pre[i]=xi_ind_pre[i]+h*function_xi_ind(xi_ind_pre[i],tau_ext,xi_ind[i]/sqrt_h,xi_all/sqrt_h, eta_c );
//xi_all_pre[i]=xi_all_pre[i]+h*function_xi_all(xi_all_pre[i],tau_ext,xi_ind[i]/sqrt_h,xi_all[0]/sqrt_h, eta_c );
i_e_ext[i]= function_i_e_ext(i_e_ext0, sigma_ext_e,eta_ind[i]);
i_i_ext[i]= function_i_i_ext(i_i_ext0, sigma_ext_i,eta_ind[i]);
}
for(int i=0;i<NN;i++)
{
if( step>(max_t_lat ))
{
current_arise_e[i]=current_arise_e_pre[i]+h*function_current_arise(current_arise_e_pre[i], tau_arise,r_e_lat[i][max_t_lat] );
current_decay_e[i]=current_decay_e_pre[i]+h*function_current_decay(current_decay_e_pre[i], tau_decay,current_arise_e[i] );
current_arise_i[i]=current_arise_i_pre[i]+h*function_current_arise(current_arise_i_pre[i], tau_arise,r_i_lat[i][max_t_lat] );
current_decay_i[i]=current_decay_i_pre[i]+h*function_current_decay(current_decay_i_pre[i], tau_decay, current_arise_i[i] );
}
else
{
current_arise_e[i]=current_arise_e_pre[i]+h*function_current_arise(current_arise_e_pre[i], tau_arise,r_e[i] );
current_decay_e[i]=current_decay_e_pre[i]+h*function_current_decay(current_decay_e_pre[i], tau_decay,current_arise_e[i] );
current_arise_i[i]=current_arise_i_pre[i]+h*function_current_arise(current_arise_i_pre[i], tau_arise, r_i[i] );
current_decay_i[i]=current_decay_i_pre[i]+h*function_current_decay(current_decay_i_pre[i], tau_decay, current_arise_i[i] );
}
current_decay_e_tem[i][0]= current_decay_e[i];
}
//**calculate current;
if(step>delay_step)
{
for(int i=0;i<NN;i++)
{
if(color[i]<1)
{
i_e[i]=i_e_pre[i]+h*delay_function_i_e(i, NN, i_e_pre[i], tau_e[i], i_e_ext[i], omega_ee, omega_ei, current_decay_e, current_decay_i[i],weight,current_decay_e_tem,t_tau);
i_i[i]=i_i_pre[i]+h*delay_function_i_i(i, NN, i_i_pre[i], tau_i[i], i_i_ext[i], omega_ie,omega_ii, current_decay_e, current_decay_i[i], weight,current_decay_e_tem,t_tau);
}
else
{
i_e[i]= i_e_ext0;
i_i[i]= i_i_ext0;
}
}
}
else
{
for(int i=0;i<NN;i++)
{
if(color[i]<1)
{
i_e[i]=i_e_pre[i]+h*function_i_e(i, NN, i_e_pre[i], tau_e[i], i_e_ext[i], omega_ee, omega_ei, current_decay_e, current_decay_i[i],weight);
i_i[i]=i_i_pre[i]+h*function_i_i(i, NN, i_i_pre[i], tau_i[i], i_i_ext[i], omega_ie, omega_ii, current_decay_e, current_decay_i[i], weight);
}
else
{
i_e[i]= i_e_ext0;
i_i[i]= i_i_ext0;
}
}
}
//**convey the latency delay data
for (int i = 0; i < NN; i++)
{
for (int j = max_t_lat; j > 0; j--)//** be careful of the j range
{
r_e_lat[i][j] = r_e_lat[i][j - 1];
r_i_lat[i][j] = r_i_lat[i][j - 1];
}
}
//** convey the delay data
for (int i = 0; i < NN; i++)
{
for (int j = max_t_tau; j > 0; j--)//** be careful of the j range
{
current_decay_e_tem[i][j] = current_decay_e_tem[i][j - 1];
}
}
// if(step==500 )
// {
// for(int i=0;i<NN;i++)
// {
// // i_e[i]=-6.7+0.2/NN*i;
// i_e[i]=-0;
// //printf("%d,%f\n",i,i_e[i]);
// i_i[i]=-1;
// }
// }
//**save data
if(step>500&step%100==0)
{
for(int i=0;i<NN;i++)
{
if(color[i]<1)
{
fprintf(fp,"%f,",i_e[i]);
fprintf(fp_i,"%f,",i_i[i]);
//fprintf(fp_xi_g,"%f,",xi_all_pre[i]);
}
}
fprintf(fp," \n" );
fprintf(fp_i," \n" );
fprintf(fp_xi_g,"%f\n",xi_all_pre);
}
/*Print progress rate*/
if(step>0&step%10000==0)
{
printf("The program has been run %.2f%%\n",(step*100.0/max_step));
endTime = clock();
std::cout << "Time already passed: " <<(double)(endTime - startTime) / CLOCKS_PER_SEC << "s" << std::endl;
}
}
fclose(fp);
fclose(fp_i);
fclose(fp_xi_g);
return 0;
}
Computing file changes ...
