% ANC mechanism
clear; clc; close all

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

fe=0.75;

J = 6;          % [0.6 - 12] 
% K = 35;              
Vr = 50;
r = 0.01;
Fi = 0.15;      % Parameter shaping the dynamics of E and Cb (s)
T = 0.2;        % Capillary transit time
Tg = J*Vr*r*T/Eo;       % Time constant of g (s) -> range (0.12 : 2.4)
Go = 0.1;       % Baseline ration of tissue and arterial plasma oxygen concent
C_Bo = Go - Eo/(log(1-Eo/(1-Go)));
Vo = 0.02;      % Resting blood volume fraction

% BOLD paramters
k1 = 7*Eo;      % Constant in the observation equation
k2 = 2;         % Constant in the observation equation
k3 = 2*Eo-0.2;  % Constant in the observation equation


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

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

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


Ks = [0.15, 0.35, 0.55];     % Neuro-CMRO2 coupling gain  (Star parameter)

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

nei=nei';
nei(nei<n0)=n0;

n = fe*nei(1,:) + (1-fe)*nei(2,:);
n = [ones(1,100),n(1:24901)];

for i=1:3   
K=Ks(i);         % Neuro-CMRO2 coupling gain  (Star parameter)

%%  HRF Zheng
HRF_Zo = [1;1;0;1;1;Eo;C_Bo;Go];
opt_hrf = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-4 1e-4 1e-4 ...
1e-4 1e-4 1e-4],'MaxStep',dt);
[~,HRF_Z] = ode15s(@(t,HRF_Z) anc_hrf_f(t, dt, HRF_Z, n, e_, Ts, Tf, Eo, Go, C_Bo, Tg, Fi, K, Tmtt, Alpha, Beta,b, Tw),...
    tspan, HRF_Zo, opt_hrf);

HRF_Z = HRF_Z';

fa = HRF_Z(1,:);    % Relative cerebral blood flow in the arteriole
v = HRF_Z(2,:);     % Relative cerebral blood volume 
s = HRF_Z(3,:);     % Vasoactive signal
q = HRF_Z(5,:);     % Relative amount deoxy-hemoglobin balloon/windkessel compart
a = HRF_Z(4,:);
E = HRF_Z(6,:);
C_B = HRF_Z(7,:);
g = HRF_Z(8,:);

%% 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*i+105)/255],"LineWidth",3);
xticklabels([])
grid on
title ('CBF')

subplot(3,1,2)
plot(tspan,v,'Color',[0 0 (50*i+105)/255],"LineWidth",3);
xticklabels([])
grid on
title ('CBV')

subplot(3,1,3)
plot(tspan,y,'Color',[0 0 (50*i+105)/255],"LineWidth",3);
hold on
xlabel("Time (s)");
grid on
title ('BOLD')

end

