sfo-ang-tnf-kca

%% Ang II + TNF Model with KCa and Ca2+ dynamics:
% Replacing the IKS with IKCa and changing internal Calcium dynamics. Calcium enters
% the cell during spiking via N-type Ca2+ channels (tied to ICa).
clear
  
tix=0;
for yu = 1
    tix=tix+1;
    
%% Single Recording Trace:
% Time:
dt = 0.01; % Step Size
t=0; % Start time (ms)
timesec = 20; % End time in s
tfinal = timesec*1000; % End time (ms)
Tstep =(0:dt:tfinal); % Time Array in ms
 
%% Initializing Values:
V=-60;       % Membrane Potential (mV)
nA=0;        % IK Activation
sA=0;        % INa Activation
kA=0;        % INa Inactivation
mNaP=0;      % INaP Activation
hNaP=0;      % INaP Inactivation
mNCa=0;      % N-type Calcium
mKA=0;       % A-type Potassium
hKA=0;       % A-type Potassium
 
%% Passive Cell Properties
% Cell Morphology:
CellSize = 10;                       % Diameter of Model Cell in microns
CellArea = 4*pi*((CellSize/2).^2);   % Area of Cell with Spherical Shape
% Cell Capacitance:
Cr = 5;                     % Real Cell Capacitance in pF
C = (Cr/CellArea)*100;      % Model Capacitance in uF/cm^2
% Cell Resistance:
Rr = 1;                     % Real Cell Resistance in Gigaohms
R = (Rr * CellArea)*10;     % Model Cell Resistance in ohms*cm^2
% Leak Conductance:
gLeak = (1/R)*1000;         % Leak Conductance to Set Input Resistance in mS/cm^2
Er = -65;                   % Reversal Potential of Leak to dictate RMP
% Proper Cell Size for measurement conversion
SizeCM = CellSize/10000;                      % cell diameter in cm
CellAreaReal = 4*pi*((SizeCM/2).^2);          % Area of Cell with Spherical Shape (cm^2)
 
%% Active Cell Properties
% Reversal Potentials (mV) from Lab:
ENa = 107;
EK = -88;
ENSCC = -35;
ECa = 120;

% Conductances (in mS/cm^2):
% Potassium 
gK = 100;                         % Conductance of DR K+                                                                                                                  ;                          % Conductance of slow act. K+-like channel
%gA = 2;                           % Conductance of A-type K+
gKCa = 0.2;                      % Conductance of KCa
% Sodium
gNa = 150;                        % Conductance of NaT
gNaP = 0.13;                      % Conductance of NaP
% Non-selective Cationic Current
%gNSCC = 0.2;                    % Conductance of Non-selective Cationic Cond.
% Calcium:
%gNCa = 0.3;                     % Conductance of N-type Calcium

% Calcium Dynamics:
F = 96.485;       % Faraday's constant
CaI = 0;          % Initiatilize Ca2+: 100 nM or 0.1 uM at rest
tauCaI = 3000;    % Rate of Ca2+ decay/clearance from neuron
SAV = (3/10);     % Surface area-to-volume ratio = s/r where s = 3 for sphere and r = 10um, divided by 1000 for cm
CaH = 0.4;        % Half activation of KCa channel (Ca dependent)
CaS = 0.1;        % Slope of activation 

% Time Constants (IK is in the loop):
tau_sA = 0.1;           % Time Constant for m gate (Sodium Act.)
tau_kA = 1;             % Time Constant for h gate (Sodium Inact.)
tau_mNaP = 5;           % Time Constant for m gate (Persistent Na+)
tau_hNaP = 50;          % Time Constant for h gate (Persistent Na+)
tau_mNCa = 5;           % Time Constant for m gate (N-type Calcium) 
tau_mKA = 5;            % Time Constant for m gate (A-type Potassium) 
tau_hKA = 30;           % Time Constant for h gate (A-type Potassium)


%% Zeros Vectors
Vm=zeros(1,length(Tstep));
Iz=zeros(1,length(Tstep));
CAI=zeros(1,length(Tstep));
ZZCA=zeros(1,length(Tstep));

%% Looping Code
% Solving DE With Eulers Method:
tidx=0;
for z=Tstep
    tidx=tidx+1;
    
%% Injected Current (1 pA = 0.3183 uA/cm^2 & 10 pA = 3.18 uA/cm^2);
I = -1;


%% Changing conductance parameters to mimic Ang II application
if (z>=2000) && (z<18000)
        gNSCC = 0.205;       % Percent Increase =  __%
        gA = 1.2;            % Percent Decrease =  __% 
        gNCa = 0.2;    
else
        gNSCC = 0.1;  
        gA = 3;
        gNCa = 0.2; 
end

%% Time Constant for IK (in ms):
tau_nA = 7.2-(6.4/(1+exp((V+28.3)/-19.2)));   % Time Constant for activation of IK

%% Current Equations:
% Potassium Currents:
    % Potassium Channel (IK):      
    nA0 = 1/(1+exp((V-2)/-8));         
    nA = nA + dt*((nA0-nA)/tau_nA);
    IK = gK*nA^4*(V-EK);
    % Transient Potassium Current (IA):
    mKA0 = 1/(1+exp((V+44)/-18));         
    mKA = mKA + dt*((mKA0-mKA)/tau_mKA);
    hKA0 = 1/(1+exp((V+61)/8));  
    hKA = hKA + dt*((hKA0-hKA)/tau_hKA);
    IA = gA*mKA^3*hKA*(V-EK);
% Sodium Currents:
    % Sodium Channel (INa):               
    sA0 = 1/(1+exp((V+31)/-6.11));  
    sA = sA + dt*((sA0-sA)/tau_sA);
    kA0 = 1/(1+exp((V+62)/6.15));
    kA = kA + dt*((kA0-kA)/tau_kA);
    INa = gNa*sA^3*kA*(V-ENa);
    % Persisent Na+ Current (INaP):  
    mNaP0 = 1/(1+exp((V+55)/-4));       
    mNaP = mNaP + dt*((mNaP0-mNaP)/tau_mNaP);  
    hNaP0 =1/(1+exp((V+45)/6.1));           
    hNaP = hNaP + dt*((hNaP0-hNaP)/tau_hNaP);
    INaP = gNaP*mNaP^3*hNaP*(V-ENa);  
% Calcium Currents:
    % N-Type Calcium Current
    mNCa0 = 1/(1+exp((V+14)/-5.8));             
    mNCa = mNCa + dt*((mNCa0-mNCa)/tau_mNCa);
    ICa = gNCa*mNCa^2*(V-ECa);
% Calcium-activated K+ Currents:
    % Internal Calcium:
    CaI = CaI + dt*(SAV*((-ICa)/F)-((CaI)/tauCaI));
%     if CaI < 0.1
%        CaI = 0.1;    % Don't drop below baseline 0.1 uM Ca2+
%     end
    % SK Channel (no voltage dependence, only Ca2+ dependent):
    zCa = 1/(1+exp(-((CaI-CaH)/CaS)));
    IKCa = gKCa*zCa*(V-EK);
% Leakage Currents:
    IL = gLeak*(V-Er);
% Non-selective Cation Current:
    INSCC = gNSCC*(V-ENSCC);
% Noise:
    IN=9*randn;
 
%% Voltage Eqn: 
    V = V + ((dt/C)*(I + IN - IL - INa - IK - INSCC - INaP - IA - ICa - IKCa));
     
% Convert Current to nA
    IKAPP = (IK+IA) * CellAreaReal;               % Total K+ Current in uA
    IKreal = IKAPP*1000;                          % Total K+ Current in nA
    INAA = (IKCa)*CellAreaReal;                   % IKCa in uA
    INAreal = INAA*1000;                          % IKCa in nA
    
    %% Vector of Values:
    Vm(tidx)= V;
    Iz(tidx) = INAreal;
    CAI(tidx)= CaI;
    ZZCA(tidx) = zCa;
       
end

%% Spikes Times & CV:
NEWNEW = Vm(1:length(Vm));
TIME = Tstep(1:length(Tstep));
FF = zeros(1,length(NEWNEW));
times = zeros(1,length(NEWNEW));
pp = 0;
for nn = 1:length(NEWNEW)
pp=pp+1;
if NEWNEW(pp) >= 15
  FF(pp) = 1;
else
  FF(pp) = 0;
end

if FF(pp)==1
times(pp)=nn; 
else
times(pp)=0;
end
end

% If vector is all 0's CV =0:
if times(:,:)<=0
    CV = 0;
else

% Calculate ISI Times:
L = times(times~=0);                % Vector of all the spike times
LL = zeros(1,length(L));
tidl = 1;
for y=2:length(L)
    tidl = tidl+1;
LL(tidl) = L(tidl) - L(tidl-1);     % Vector of interspike intervals (has double points though)
end

tidr = 0;
for rr=1:length(L)
tidr = tidr + 1;

if LL(tidr) <= 1       % If ISI is less than 1 (two points on the same spike -- both above 15 mV)
    L(tidr) = 0;
else
    L(tidr) = L(tidr);
end
    
woohoo = L(L~=0);      % Woohoo is the proper spike times (excluding any double points from the same AP).
end

%% Calculate CV
time=2:length(woohoo);
ISIT=zeros(1,length(time));
N=zeros(1,length(time));

n0=1;
tidt=0;
for p=time
    tidt=tidt+1;
    
    ISI=woohoo(n0+1)-woohoo(n0);
    n0=n0+1;
    
    ISIT(tidt)=ISI;
    N(tidt)=n0;
end

ISIArray=transpose(ISIT);

% CV = STD of ISIs/ Mean ISIs :
CV = std(ISIArray)/mean(ISIArray);

end

%% End of Running over many traces:
% CVAVG(tidf) = CV;
% end

% CVAvg = mean(CVAVG)
% STD = std(CVAVG)

%% Plotting MP as a Histogram:
% VMM=Vm';
% figure
% H = histogram(VMM,200);

%% Spike Train & Firing Frequency: 
% Vector of 0&1's at proper times:
SPIKEPLACE = zeros(1,length(Tstep));
sptimes = woohoo;
sptimes(1,length(woohoo)+1)=0;
hhh=0;
hfh=0;
for pvect = 1:length(Tstep)
   hhh = hhh+1;
    
if pvect == sptimes(hfh+1)   
   SPIKEPLACE(hhh)=1;         % SPIKEPLACE contains a 1 when there is a spike
   hfh = hfh+1;
else
   SPIKEPLACE(hhh)=0;
   hfh = hfh;
end
end

% Calculate FF in Hz (how many spikes/second)
seg = 100000;    % Every 1s 
tiff = 100000;
tif = 1;
tifd = 0;
for tt = 1:seg:length(Tstep)
    tifd = tifd+1;   
    firefreq = sum(SPIKEPLACE(tif:tiff));  % Sum the array vector every 1s
    tif = tif+seg;
    tiff = tiff+seg;
    FiringFreq(tifd)=firefreq;  % Firing frequency of the neuron every 1 second (if 60s recording, 60 FF points)
    
    if tif == length(Tstep)
    break    
    end
end

MeanFiringFrequency = mean(FiringFreq);   % Mean FF over entire trace.

%% Calculate CV and FF:
% CALCULATECV(tixf)= CV;
% CALCULATEFF(tixf)= MeanFiringFrequency;
% VHALFF(tixf) = vhalf;

MeanCaI = mean(CAI(200000:1800000));
CAIFORPRISM(tix,:) = CAI;

LARGECV(tix)=CV;
LARGEFF(tix)=MeanFiringFrequency;
LARGECAI(tix)=MeanCaI;
end

MeanLargeCV = mean(LARGECV);
MeanLargeFF = mean(LARGEFF);
STDFF = std(LARGEFF)/sqrt(length(LARGEFF));
% LARGERFF(tiw,:)=LARGEFF;

% hold on
% data_mean=mean(CAIFORPRISM,1);
% plot(data_mean)



%% Graphs: 
CVV = num2str(CV);
FFF = num2str(MeanFiringFrequency);
%CALCIUMCON = num2str(MeanCaI);
% % 
figure
sh(1)=subplot(3,1,3);
plot(Tstep,Vm,'k-')
xlabel('Time (ms)')
ylabel('Voltage')
% % % ylim([-80 80])
% IKCa Graph:
sh(1)=subplot(3,1,2);
plot(Tstep,Iz,'k-')
xlabel('Time (ms)')
ylabel('IKCa')
% Calcium Graph:
sh(1)=subplot(3,1,1);
plot(Tstep,CAI,'k-')
xlabel('Time (ms)')
ylabel('[Calcium]')
linkaxes(sh,'x')
title(['CV=',CVV,' & Mean FF=',FFF,'Hz'])