% ECI mechanism
clear; clc;

%% Parameters
e_ = 0.28;      % Neurovascular coupling gain  -> range (0.2 : 0.5)
Ts = 1.1;       % Time constant of vasoactive signal (s) range (0.9: 1.2)
Tf = 1.6;       % Autoragulation time constant (s) range (1.4 : 2.2)
Tmtt = 1.8;     % Total tissue mean tansit time (s) -> range (1.2 : 2.4)
Beta = 2;       % Diminished reserved volume exponent
Alpha = 2;      % Exponent of Poiseville's law 
Eo = 0.5;       % Steady state oxygen extraction fraction
Go = 0.1;       % Steady state value of g
Tw = 6.68;      % Time constant of a (s)
b = 7.6;        % Gain modeling the influence of the variation of CBF a (s)
n0 = 0.4;          % minimun allowed of n(s);
T_e = 10e-3;    % Neuronal exitation time constant region 1
% T_i = 30e-3;       % Neuronal inhibitory time constant region 1
ce = 1;        % Excitatory-Inhibitory conectivity region 1                
ci = 2;       % Inhibitory-Excitatory conectivity region 1 

ceu = 1;        % Input strength to excitatory state in region 1
ciu = 2;      % Input strength to inhibitory state in region 1


fe=0.75;


Vo = 0.02;      % Resting blood volume fraction
v0 = 40.3;      % Ratio of intra-vascular to extra-vascular fMRI signal 
                % contributions (Hz)
Te = 45e-3;     % Echo time (s)
Epsilon = 1.32; % Frequency offset at the surface of a blood vessel for 
                % fully deoxygenated blood for Gradient Echo (GE)
                % sequence at 1.5T (Hz)
r0 = 15;        % Sensitivity of changes in intra-vascular signal relaxation
                % rate with changes in oxygen saturation for GE sequence
                % at 1.5T                

k1 = 4.3*v0*Eo*Te;          % Constant in the observation equation
k2 = Epsilon*r0*Eo*Te;      % Constant in the observation equation
k3 = 1-Epsilon;             % Constant in the observation equation


% Temporal parameters
maxt = 25;
dt = 1e-3;
tspan = 0:dt:maxt;

%% Input spikes
u0 = 28;
input = IED_function(tspan,u0);


Ti_all = [30e-3, 100e-3, 200e-3];       % Star parameter

for j=1:3 
%% Simulating neuronal activity
T_i = Ti_all(j);                        % Star parameter

n0_ei = [1;1];
opt_neuro = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4],'MaxStep',dt);
[~,nei] = ode15s(@(t,nei) eci_neuronal(t,dt,input,nei,n0,T_e,T_i,ce,ci,ceu,ciu),...
    tspan, n0_ei, opt_neuro);

nei=nei';
n = fe*nei(1,:) + (1-fe)*nei(2,:);

%% HRF
HRFo = [1;1;1;0;1;1];
hrf_opt = odeset('MaxStep',dt);
[~,HRF] = ode15s(@(t,HRF) eci_hrf_f(t, dt, HRF, n,e_, Ts, Tf, Eo, Go, b, Tw, Tmtt, Alpha, Beta), tspan, HRFo, hrf_opt); 

HRF = HRF';

fa = HRF(1,:);    % Relative cerebral blood flow in the arteriole
r = HRF(2,:);     % Relative arteriole resistance
v = HRF(3,:);     % Relative cerebral blood volume 
s = HRF(4,:);   % Vasoactive signal
a = HRF(5,:);    % Relative pressure to volume relation constant
q = HRF(6,:);   % Relative amount deoxy-hemoglobin balloon/windkessel compart


%% Calculating BOLD Response

   y = Vo*(k1*(1-q) + k2*(1-(q./v)) + k3*(1-v)); 
   

%% ploting section

    subplot (3,1,1)  
    plot (tspan,fa, 'Color',[0 0 (50*j+105)/255],'LineWidth',3);
    xticklabels([ ])
    grid on;
    hold on
    
    subplot (3,1,2)
    plot (tspan,v, 'Color',[0 0 (50*j+105)/255],'LineWidth',3);
    xticklabels([ ])
    grid on;
    hold on

    subplot(3,1,3)
    if j==1
        plot (tspan,y(1,:),'Color',[0 0 (50*j+105)/255],'LineWidth',3);
        xlabel("Time (s)");
        grid on;
        hold on
    end
    if j==2
        plot (tspan,y(1,:),'Color',[0 0 (50*j+105)/255],'LineWidth',3);
        xlabel("Time (s)");
        grid on;
        hold on
    end
    if j==3
        plot (tspan,y(1,:),'Color',[0 0 (50*j+105)/255],'LineWidth',3);
        xlabel("Time (s)");
        grid on;
        hold on
    end
    
end










