%
% Ref. 
% G. Valenza et.al. "A neuron–astrocyte transistor-like model for
% neuromorphic dressed neurons", Neural Networks, Volume 24, Issue 7, September 2011, Pages 679-685
%
%

clear all
clc
warning off,
% define constant parameters

k1=78.59;       k2=136.07;           c1=0.14768;       c2=0.022674;             n1=-2.0521;       n2=-0.86603;

a00=20.791;   a01=-105.7;   a02=257.04;   a03=-248.63;   a04=82.978;
a10=-0.20103;   a11=1.9447;   a12=-3.5931;   a13=2.0715;
a20=-2115.3;   a21=16048;   a22=-3054.5;   a23=257.21;   a24=-8.5669;
a30=2531.3;   a31=2502.6;   a32=-148.85;
a40=1.3381;   a41=-0.28508;   a42=0.018957;
a50=3.4867;   a51=33.483;   a52=-51.257;   a53=26.083; % a0 di Ifin=a0+a1In  
a60=0.00050985;   a61=0.0018472;      % a1 di Ifin=a0+a1In
a70=45736;   a71=-2853;    %delay 2
a80=7.2386;   a81=-0.51501; %overshoot
a90=0.0069871;   a91=0.0083831;   %f
a100=1460.8;  a101=-234.4;  a102=9.4377;   %tau


% def function


r=input('Insert the rIP3 value: ');
In=input('Insert the pre-synaptic current value In [pA]: ');

res=0.001;
tfin=input('Insert the simulation time [s]: ');
t1=(0:res:tfin); %time in sec

s1=k1*((r-c1)^(1/n1));

s2=k2*((r-c2)^(1/n2));

%%%%%%%%%%%%%%

if In<=s1
        Iastro1=zeros(1,length(t1));
end


if (In>s1) && (In<=s2)
    
    h=[];
    k3=a00 + (a01*r) + (a02*(r^2)) + (a03*(r^3)) + (a04*(r^4));

    k4=a10 + (a11*r) + (a12*(r^2)) + (a13*(r^3));

    A1=(k3*((In-s1)^k4));           % /In;

    d1=a20 + (a21*A1) + (a22*(A1^2)) + (a23*(A1^3)) + (a24*(A1^4)); % delay in msec
    d1=floor(d1);

    T=a30 + (a31*A1) + (a32*(A1^2)); %in msec
    dc=a40 + (a41*A1) + (a42*(A1^2));

    m=pi/(dc*T);

    T1=dc*T;
    T1=round(T1);
    T2=T-T1;
    T2=round(T2);
    
    while(length(h)<length(t1))
        h=cat(2,h,[m*t1(1:T1)*1000,zeros(1,T2)]);
    end
    Iastro1=[zeros(1,d1),A1*sin(h(1:(length(t1)+1-d1)))];
       Iastro1(end)=[];
 
end

 


if In>s2
    
    k5=a50 + (a51*r) + (a52*(r^2)) + (a53*(r^3));

    k6=a60 + (a61*r);

    Ifin=k5 + (k6*In);

    d2=a70 + (a71*Ifin); % delay in msec
    d2=floor(d2);
    A2=a80 + (a81*Ifin);

    f=a90 + (a91*Ifin);

    tau=a100 + (a101*Ifin) + (a102*(Ifin^2));
    
    Iastro1=[zeros(1,d2),Ifin+(A2*exp(-t1(1:end+1-d2)./tau).*sin(2*pi*f*t1(1:end+1-d2)))];
             Iastro1(end)=[];
 
 end
  
figure
plot(t1,Iastro1),xlabel('Time [s]'),ylabel('Iastro [pA]')

