% ABS mechanism
clear; clc;

%% Paramentes
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;10e-3];    % Neuronal exitation time constant region 1 & 2 
T_i = [30e-3;30e-3];       % Neuronal inhibitory time constant region 1 & 2
ce = [1;1];        % Excitatory-Inhibitory conectivity region 1 & 2                
ci =  [2;2];       % Inhibitory-Excitatory conectivity region 1 & 2
ceu = [1;0];        % Input strength to excitatory state in region 1 & 2
ciu = [0;0];        % Input strength to inhibitory state in region 1

fe=0.75;
K = 0.05;       % Neuro-CMRO2 coupling gain                
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)

% BOLD parameters
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

tauI = 1e-3;        % sinductance time constant

Ra0 = 0.39;         % baseline arteriole resistance of each region

% Mass matrix
Mass = blkdiag(ones(2),eye(10));

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


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

% Square pulse
% input = zeros(size(tspan));
% input(2:1001) = ones(1,1000);

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

nei=nei';

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

for ii = 1:3
%% HRF
RI = 0.13+(ii-1)*0.07;      % resistance of the feeding arteriole (Star parameter)
Rw0 = 1 - (2*RI + Ra0);     % windkessel resistance of each region

HRFo = [1;1;1;1;1;1;0;0;1;1;1;1];

opt_hrf = odeset('Mass',Mass,'RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-4 1e-4 1e-4 ...
    1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4],'MaxStep',dt);

[t,HRF] = ode15s(@(t,HRF) abs_sim_f(t, dt, n, HRF, e_, Ts, Tf, Eo, Go, b, Tw, Tmtt, Alpha, Beta, tauI, Ra0, Rw0, RI),...
    tspan,HRFo,opt_hrf); 

HRF = HRF';

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

subplot(3,1,2)
plot(t,v(1,:),'Color',[(50*ii+105)/255 0 0],"LineWidth",3)
hold on
plot(t,v(2,:),'Color',[0 0 (50*ii+105)/255],"LineWidth",3)
hold on
grid on
title('CBV')

subplot(3,1,3)
plot(t,y(1,:),'Color',[(50*ii+105)/255 0 0],"LineWidth",3);
hold on
plot(t,y(2,:),'Color',[0 0 (50*ii+105)/255],"LineWidth",3);
hold on
grid on
title('BOLD')
end





