% Generates multiple data files of 55 hours of stimulation of HHSIP model
clear all;
% close all;
clc
% This simulation does not require any sub-functions
% 1D slow inactivaion (gM=0)+stochasticity
% HH rates and Capacitance adjusted, Also s rates (gamma and delta)
load('HHSIP_params_7.5_7.7_7.9_8.1_8.3_dt_5e-4.mat');
intermittent=1; %if intermittent mode is required, set as 1, if transient mode required, set as 0
I_array=[7.5 7.7 7.9 8.1 8.3];% [microamper]
D=0; % current diffusio constant [(microA)^2/milisec]
f_array=[10 20 30 40 50];% [Hz]
L_I=length(I_array);
T_array=1e3./f_array; %[ms]
L_f=length(f_array);
Time=55*1e3*3600; %[msec] length of simulation
% cell_AP=cell(L_I,L_f);cell_s1=cell_AP;cell_s2=cell_AP;cell_I=cell_AP;
% cell_Latency=cell_AP;cell_Amplitude=cell_AP; %sampling arrays
% cell_stim_time=cell_AP;
%% HH parameters
N=1e6; %noise level
phi_HH=2; % make model faster
VNa=50; VK=-77;VL=-54; %[mV]
gNa=120;gK=36;gL=0.3;gM=0.05*gK; %[mS]
Cm=1/phi_HH; %[microFarad] - not 1 Cm as in HH, since spike time is shorter
dt=0.01/phi_HH; %[msec] time step
t_0=1/phi_HH; %[ms] pulse width
I_step=sqrt(D*dt);
sqrt_dt_N=sqrt(dt/N);
%% Slow kinetics parameters
phi_s=1/20; %s is slower, to get relaxation as in Gal 2010, and
amp=0.0034*3;sigma=1/0.3;Vhalf=-17; % makes inactivation faster to make
sigma2=15;
%% Main Simulation
cell_AP=cell(L_I,L_f);cell_s1=cell_AP;cell_s2=cell_AP;
tic
for jj=5%1:L_I
for ii=1%1:L_f
f_in=f_array(ii);
T=1e3/f_in; %[msec]
stim_num=round(Time/T);
stim_time=T*(0:stim_num-1);
AP=zeros(stim_num,1);s1=AP;s2=AP;Latency=NaN*AP;Amplitude=NaN*AP;I_sample=NaN*AP; %sampling arrays
RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(100*clock))); % so rand won't generate the same pattern each time...
defaultStream = RandStream.getDefaultStream;
%set initial conditions and intialize arrays
% if intermittent==1
% s_initial=params(jj,2)+0.001; %start just above threshold for intermittent mode
% else
% s_initial=1;
% end
I0=I_array(jj);
% y0=[ -66.4379 0.0446 0.0446 0.0446 0.2959 0.2959 0.2959 0.2959 0.6451 s_initial 0.77];%initial condition when s is approximatly at steady state
y0=[ -66.4379 0.0446 0.0446 0.0446 0.2959 0.2959 0.2959 0.2959 0.6451 0.89 0.85];%initial condition when s is approximatly at steady state
y=y0; %intial conditions
dy2dt=zeros(1,11); dy=dy2dt;
alpha=zeros(1,5); %recovery rates
beta=alpha; %inactivation rates
V_threshold=-10; %[mV] AP is defined if V crosses V_threshold upwards
stim_flag=0;
delay=0;
Vmax=-inf;
cycle=round(T/dt); %note that this is an approximation of the resulted T, not the real one
Pulse_width=round(t_0/dt);
for pp=1:stim_num
I_sample(pp)=I0;
s1(pp)=y(10);
s2(pp)=y(11);
delay=-dt;
stim_flag=1;
for kk=1:cycle
if stim_flag==1 %sample
delay=delay+dt;
if (y(1)>V_threshold)
AP(pp)=1;
Vmax=y(1);
stim_flag=0;
end
end
if Vmax>-200
if (y(1)<Vmax) %sample
Latency(pp)=delay;
Amplitude(pp)=Vmax;
Vmax=-inf;
else
Vmax=y(1);
delay=delay+dt;
end
end
I=I0*(kk<Pulse_width);
%% SO_hhx_Langevin inserted here
V=y(1);
beta(1)= phi_HH*0.125.*exp(-(V+65)./80); %n
beta(2)= phi_HH*4.*exp(-(V+65)./18); %m
beta(3)= phi_HH*1/(exp(-0.1*(V+35))+1); %h
beta(4)= phi_s*amp./(exp(-(V-Vhalf)/sigma)+1); % s1 - different then Fleidervish1996 !!
s2_inf=1./(1+exp(-(V+35)/10));
tao_s2=1e3./(3.3*exp((V+35)/sigma2)+exp(-(V+35)/20));
beta(5)=(1-s2_inf)./tao_s2; %s2- orignally from Yamada1989
alpha(1)=phi_HH*0.01*(V+55)./(1-exp(-0.1*(V+55)));%n
alpha(2)=phi_HH*0.1*(V+40)./(1-exp(-0.1.*(V+40)));%m
alpha(3)=phi_HH*0.07.*exp(-(V+65)./20); %h
alpha(4)=phi_s*0.001.*exp(-(V+85)./30); %s1
alpha(5) = s2_inf./tao_s2; %0.0034./(exp(-0.1*(V+17))+1);- orignally from Yamada1989
dy2dt(1)=(gNa.*(y(2).*y(3).*y(4)).*y(9).*y(10).*(VNa-y(1))+(y(5).*y(6).*y(7).*y(8)).*(gK+gM.*y(11)).*(VK-y(1))+gL.*(VL-y(1))+I)./Cm;% dVm/dt
dy2dt(2)=(1-y(2)).*alpha(2)-y(2).*beta(2); % dm/dt
dy2dt(3)=(1-y(3)).*alpha(2)-y(3).*beta(2); % dm/dt
dy2dt(4)=(1-y(4)).*alpha(2)-y(4).*beta(2); % dm/dt
dy2dt(5)=(1-y(5)).*alpha(1)-y(5).*beta(1); % dn/dt
dy2dt(6)=(1-y(6)).*alpha(1)-y(6).*beta(1); % dn/dt
dy2dt(7)=(1-y(7)).*alpha(1)-y(7).*beta(1); % dn/dt
dy2dt(8)=(1-y(8)).*alpha(1)-y(8).*beta(1); % dn/dt
dy2dt(9)=(1-y(9)).*alpha(3)-y(9).*beta(3); % dh/dt
dy2dt(10)=(1-y(10)).*alpha(4)-y(10).*beta(4); % ds1/dt
dy2dt(11)=(1-y(11)).*beta(5)-y(11).*alpha(5); % ds2/dt - notice the switch in alpha and beta!
noise(2)=sqrt(((1-y(2)).*alpha(2)+y(2).*beta(2))); % dm/dt
noise(3)=sqrt(((1-y(3)).*alpha(2)+y(3).*beta(2))); % dm/dt
noise(4)=sqrt(((1-y(4)).*alpha(2)+y(4).*beta(2))); % dm/dt
noise(5)=sqrt(((1-y(5)).*alpha(1)+y(5).*beta(1))); % dn/dt
noise(6)=sqrt(((1-y(6)).*alpha(1)+y(6).*beta(1))); % dn/dt
noise(7)=sqrt(((1-y(7)).*alpha(1)+y(7).*beta(1))); % dn/dt
noise(8)=sqrt(((1-y(8)).*alpha(1)+y(8).*beta(1))); % dn/dt
noise(9)=sqrt(((1-y(9)).*alpha(3)+y(9).*beta(3))); % dh/dt
noise(10)=sqrt((((1-y(10)).*alpha(4)+y(10).*beta(4)))); % ds1/dt - notice phi!
noise(11)=sqrt(((1-y(11)).*beta(5)+y(11).*alpha(5))); % ds2/dt - notice the switch in alpha and beta!
I0=I0+I_step*(rand-0.5);
if I0>8.3
I0=8.3;
end
y=y+ dy2dt*dt+randn(1,11).*noise.*sqrt_dt_N;
%% here SO_hhx_Langevin ends
end
end
cell_AP(jj,ii) ={AP};cell_s1(jj,ii) ={s1};cell_s2(jj,ii) ={s2};
clear AP s1 s2 stim_time
% save(['HHSIP_N=1e6_f=10-50_I_7.5-8.3_T_1e7.mat' ]);
save(['HHSIP_N=1e6_f=10_I_8.3_Time_1e6.mat' ]);
end
end
toc