% NDA 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
Tw = 6.68;      % Time constant of a (s)
b = 7.6;        % Gain modeling the influence of the variation of CBF a (s)
n0 = 0.3;       % minimun allowed of n(s);

% T_e = [10e-3;3];        % Neuronal exitation time constant region 1 & 2 
                        % for NDA -> T_e2 range (1 : 10)
T_i = [30e-3;30e-3];    % Neuronal inhibitory time constant region 1 & 2

ce = [1;0.01];        % Excitatory-Inhibitory conectivity region 1 & 2 for NDA 
ci = [2;2];        % Inhibitory-Excitatory conectivity region 1 & 2

ceu = [1;1];        % Input strength to excitatory state in region 1 & 2
ciu = [0;2];        % Input strength to inhibitory state in region 1 & 2
                    % for NDA and remote-LCI
                    
fe = [0.75,1];        % fraction of inhibitory neurons

Fi = 0.15;      % Parameter shaping the dynamics of E and Cb (s)
Go = 0.1;       % Steady state value of g
Tg = 1.6;       % Time constant of g (s) -> range (0.12 : 2.4)
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 = 320;
input = IED_function(tspan,u0);

% input = zeros(1,length(tspan));
% input(2:1001) = ones(1,1000);


Te_all = [2 , 4, 8];    % Star parameter

for i = 1:3
    
%% Simulating neuronal activity
T_e = [10e-3,Te_all(i)];
n0_ei = [1;1;1;1];  
opt_neuro = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-4 1e-4],'MaxStep',dt);
[~,nei] = ode15s(@(t,nei) nda_neuronal(t,dt,input,nei,n0,T_e,T_i,ce,ci,ceu,ciu),...
    tspan, n0_ei, opt_neuro);

nei=nei';

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

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

HRF = HRF';

fa = HRF(1:2,:);    % Relative cerebral blood flow in the arteriole
v = HRF(3:4,:);     % Relative cerebral blood volume 
s = HRF(5:6,:);     % Vasoactive signal
a = HRF(7:8,:);    % Relative pressure to volume relation constant
q = HRF(9:10,:);   % 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(1,:), 'Color',[(20*i+195)/255 0 0],'LineWidth',3);
hold on
plot (tspan,fa(2,:), 'Color',[0 0 (20*i+195)/255],'LineWidth',3);
hold on
xticklabels([ ])
title ('CBF')
grid on;

subplot (3,1,2)
plot (tspan,v(1,:), 'Color',[(20*i+195)/255 0 0],'LineWidth',3);
hold on
plot (tspan,v(2,:), 'Color',[0 0 (20*i+195)/255],'LineWidth',3);
xticklabels([ ])
grid on;
title ('CBV')
hold on

subplot(3,1,3)
plot(tspan,y(1,:),'Color',[(20*i+195)/255 0 0],'LineWidth',3);
hold on
plot(tspan,y(2,:),'Color',[0 0 (20*i+195)/255],'LineWidth',3);
grid on
title ('BOLD')

end
