function A2ONCB_FHN
%%% In order to understand the slow platau-like oscillation arising in rd
%%% AII after LP application, we hypothesized that there is a slow Ca2+
%%% mediated oscillation in the coupled ONCB. We model the Ca oscillation
%%% as a relaxation oscillation based on FitzHugh-Nagumo model.


clc
clear all

global epsilon A vz0 vz1 vz2 vz3 vz4 w0 gamma v0 exp_amp exp_thresh exp_steep tau_w
global delta_tau v_tau v_xi delta_tauM v_tauM v_xiM

global gjONCB ie ie_tonic ie_pulse t_min t_max

global numAII armsection R capacitance            
global gbar_na gbar_M gbar_A G_gap Gpas epas ena ek      
global vhalfh_na vhalfm_na kh_na km_na vhalfm_M km_M vhalfh_A vhalfm_A vhalfw_A kh_A km_A kw_A f
global mtau_na htau_na mtau_M mtau_A 
global tau capacitanceONCB 
global t_switch t_width gbar_M_IS_1 gbar_M_IS_2 SA_hand


%%%%%%%% -------------------- Switches ---------------- %%%%%%%%%

numAII = 1;
ONCB = 10;%10;
ie = 0;
% current injection during the duration t_min to t_max
ie_tonic=0;
ie_pulse=0;%1e-8;
t_min=254;
t_max=263;

%% 
%%%%%%%% --------- FitzHugh-Nagumo parameters in ONCB --------- %%%%%%%%%

epsilon = 0.000002;
A = -0.01;
vz0 = -58;
vz1= -48;
vz2 = -40;%-30;%-17;%-27;%-30;
vz3 = -35;%-15;%-25;
vz4 = -61;%-63;
w0 = 0;
gamma = 30;%50;
v0 = -58.5;%-58;%-55;
exp_amp=0.3;
exp_thresh=-56;
exp_steep=0.2;
tau = 8000;%2000;%8000;%16000;%5000;
ONCB_area_factor=0.5;%0.64;%0.2;
time_end = 3000;
tau_w=2.5;
delta_tau=-.95;%-0.95;
v_tau=-55;
v_xi=-1;%-0.1;
delta_tauM=0;
v_tauM=-45;
v_xiM=-0.1;

v_shift=-1;
vz0=vz0+v_shift; vz1=vz1+v_shift;
vz2=vz2+v_shift; vz3=vz3+v_shift;
vz4=vz4+v_shift; 
v_tau=v_tau+v_shift; v0=v0+v_shift; exp_thresh=exp_thresh+v_shift;

%%
%%%%%%%% -------------------- General ---------------- %%%%%%%%%

ena = 50; %(mV)
ek = -77; %(mV)
Cm = 1e-3; %(mf/cm^2)

%%
%%%%%%%% -------------------- ONCB ------------------- %%%%%%%%%
v_init_ONCB = -60;
w_init = 0;

SA_ONCB= 439.8*1e-8; %(cm^2)
SA_ONCB= SA_ONCB*ONCB_area_factor; %(cm^2)

capacitanceONCB = Cm*SA_ONCB;


RmONCB = 12000;
epasONCB = -50;
gpasONCB=1/RmONCB;

if ONCB == 5
    gjONCB = 250*1e-12; %(S)
elseif ONCB == 7
%     gjONCB=  1250*1e-12; %(S)
    gjONCB=  750*1e-12; %(S)
elseif ONCB == 10
%     gjONCB=  1250*1e-12; %(S)
    gjONCB=  500*1e-12; %(S)
elseif ONCB == 0
    gjONCB = 0;
end

%%
%%%%%%%% -------------------- AII -------------------- %%%%%%%%%

v_init = -60;%-50;

armsection = 1; %11;

somadiam=25*1e-4;  %(cm)
somalength=25*1e-4;  %(cm)
armdiam=0.3*1e-4;  %(cm)
armlength=32*1e-4;  %(cm)
handdiam=2*1e-4;  %(cm)
handlength=2*1e-4;  %(cm)

SA_soma= 1963.5*1e-8; %(cm^2)
SA_arm=30.1593*1e-8; %(cm^2)
SA_hand=12.5664*1e-8; %(cm^2)

global_ra = 150; %(ohm-cm)

vhalfh_na = -49.5; %(mV)
vhalfm_na = -48; %(mV)
kh_na = 2; %(mV)
km_na = 5; %(mV)
mtau_na = 0.01; %(ms)
htau_na = 0.5; %(ms)

vhalfm_M = -40; %(mV) 
km_M = 4; %(mV)
mtau_M = 50;%(ms)

vhalfh_A = -40.5; %(mV)
vhalfm_A = -10; %(mV)
vhalfw_A= -45; %(mV)
kh_A = 2; %(mV)
km_A = 7; %(mV)
kw_A = 15; %(mV)
f=0.83;
mtau_A = 1; %(ms)

gbar_na_soma = 0;
gbar_na_arm = 0;
gbar_na_IS = 0.2; %(mho/cm^2) 

gbar_M_soma = 0;
gbar_M_arm = 0;
gbar_M_IS = 0.03;%0.035;%0.026;%0.03;%0.03;%0.015;%0.01;%*4/3; %(mho/cm2)
% gbar_M_IS = 0.03;


% variable g_m in IS
t_switch=1000;
t_width=1500;
gbar_M_IS_1=0.035;
gbar_M_IS_2=0.015;

figure(100)
nVAII=5;
for iVAII=1:nVAII
    VAII=-60+iVAII*5;
    VCB=-65:0.1:-15;
    gjONCB_plot=gjONCB;
    Vnull=  (tau*epsilon/capacitanceONCB) * gjONCB_plot*(VAII-VCB)...
        +A*(VCB-vz0).*(VCB-vz1).*(VCB-vz2)./(-VCB+vz3)./(VCB-vz4) - w0+exp_amp./(1+exp(-(VCB-exp_thresh)/exp_steep)); % v_ONCB
    wnull=(VCB  - v0)/gamma; %w
    subplot(1,nVAII,iVAII)
    plot(VCB,Vnull);
    hold all;
    plot(VCB,wnull);
    xlabel(' ONCB Voltage')
    ylabel(' w')
    title([' Nullclines V_{AII}=',num2str(VAII),' (mV)']);
    axis([-65 -30 -0.05 0.4])
    hold off
end
drawnow

% rd1 model has reduced A current
%gbar_A_soma = 0.004* 3/4;  %(mho/cm2) 
%gbar_A_arm = 0;
%gbar_A_IS = 0.08* 3/4;  %(mho/cm2) 

% wt has full A current
A_current_factor=1;
gbar_A_soma = 0.004*A_current_factor;  %(mho/cm2) 
gbar_A_arm = 0;
gbar_A_IS = 0.08*A_current_factor;  %(mho/cm2) 



% if numAII==1
%     Rm = 10000; %(ohm*cm^2)
%     epas = -40; %(mV)
% else
%     Rm = 40000; %(ohm*cm^2)
%     epas = -10; %(mV)
% end

%%%%%  Rd case, probably only for a network of AIIs
Rm = 40000; %(ohm*cm^2)
epas = -60;%-30; %(mV) % -60; %(mV)
gpas=1/Rm; %(S/cm^2)


gjCondR=  6000*1e-12;  %(S) : Mark's newer version (network oscillation) has 6000 pS or 4000pS, while the older version has 400 pS
gjCondIS= 1e-5*1e-12; %(S)
gjCondDC= 1e-5*1e-12; %(S)


R=zeros(armsection+1,1);
for num=1:armsection+1
    if num==1
        R(num)=(global_ra/(pi*(somadiam/2)^2))*somalength/2 + (global_ra/(pi*(armdiam/2)^2))*armlength/armsection/2;
    elseif num==armsection+1
        R(num)=(global_ra/(pi*(handdiam/2)^2))*handlength/2 + (global_ra/(pi*(armdiam/2)^2))*armlength/armsection/2;
    else
        R(num)=(global_ra/(pi*(armdiam/2)^2))*armlength/armsection;
    end
end


%%
gbar_na=zeros(numAII,armsection+2);
gbar_M=zeros(numAII,armsection+2);
gbar_A=zeros(numAII,armsection+2);
capacitance=zeros(numAII,armsection+2);
G_gap=zeros(numAII,armsection+2);
Gpas=zeros(numAII,armsection+2);

for ii=1:numAII
    
    for jj=1:armsection+2
        if jj==1
            gbar_na(ii,jj) = gbar_na_soma * SA_soma;
            gbar_M(ii,jj) = gbar_M_soma * SA_soma;
            gbar_A(ii,jj)= gbar_A_soma * SA_soma;
            capacitance(ii,jj) = Cm*SA_soma;
            G_gap(ii,jj) = gjCondR;
            Gpas(ii,jj) = gpas * SA_soma;
            
            
        elseif jj==armsection+2
            gbar_na(ii,jj) = gbar_na_IS * SA_hand;
            gbar_M(ii,jj) = gbar_M_IS * SA_hand;
            gbar_A(ii,jj) = gbar_A_IS * SA_hand;
            capacitance(ii,jj) = Cm*SA_hand;
            G_gap(ii,jj) = gjCondIS;
            Gpas(ii,jj) = gpas * SA_hand;

        else
            gbar_na(ii,jj) = gbar_na_arm * SA_arm/armsection;
            gbar_M(ii,jj) = gbar_M_arm * SA_arm/armsection;
            gbar_A(ii,jj) = gbar_A_arm * SA_arm/armsection;
            capacitance(ii,jj) = Cm*SA_arm/armsection;
            G_gap(ii,jj) = gjCondDC;
            Gpas(ii,jj) = gpas * SA_arm/armsection;
            
        end
        
    end
end


%%%%%% ------------------- Initial conditions --------------- %%%%%%%%%%
minit_na = 1/(1+exp(-(v_init-vhalfm_na)/km_na));
hinit_na = 1/(1+exp((v_init-vhalfh_na)/kh_na));
minit_M = 1/(1 + exp(-(v_init-vhalfm_M)/km_M));
minit_A = 1/(1+exp(-(v_init-vhalfm_A)/km_A));
hinit_A = f*(1/(1 + exp((v_init-vhalfh_A)/kh_A)))+(1-f);

initial = zeros(numAII, armsection+2, 9);
for ii = 1:numAII
        for jj = 1:armsection+2
            initial(ii,jj,1) = minit_na;
            initial(ii,jj,2) = hinit_na;
            initial(ii,jj,3) = minit_M;
            initial(ii,jj,4) = minit_A;
            initial(ii,jj,5) = hinit_A;
            initial(ii,jj,6) = hinit_A;
            initial(ii,jj,7) = v_init;
            initial(ii,jj,8) = v_init_ONCB;
            initial(ii,jj,9) = w_init;
            
        end
end



initial = reshape(initial,numAII*(armsection+2)*9,1);

time_start = 0;
%time_end = 800;
options = odeset('RelTol', 1e-8);
%options = odeset('RelTol', 1e-14);
[T,Y1] = ode15s(@solve_RelaxOscl_2, [time_start time_end],initial,options);
Y = zeros(numAII,armsection+2,9,length(T));
for i = 1:length(T)
    Y(1:numAII,1:armsection+2,1:9,i) = reshape(Y1(i,:),numAII,armsection+2,9);
end


plotcell = 1;
plotsection = 1; % 1 for soma, 3 for IS
v= Y(plotcell,plotsection,7,:);
v_fix=reshape(v,size(T));

plotsection = 3; % 1 for soma, 3 for IS
m_time= Y(plotcell,plotsection,3,:);
m_fix=reshape(m_time,size(T));
vs_time= Y(plotcell,plotsection,7,:);
vs_fix=reshape(vs_time,size(T));

plotsection = 1; % 1 for soma, 3 for IS
vCB= Y(plotcell,plotsection,8,:);
vCB_fix=reshape(vCB,size(T));

w= Y(plotcell,plotsection,9,:);
w_fix=reshape(w,size(T));

figure(10)
gbar_MIS_T=gbar_M_IS_1+(gbar_M_IS_2-gbar_M_IS_1)*.5*(1+tanh((T-t_switch)/t_width));
subplot(5,1,1), plot(T,v_fix), legend('V_{AII}')
%subplot(1,4,2), plot(T,vs_fix), legend('V_{AII} IS')
subplot(5,1,2), plot(T,m_fix), legend('m_{AII}')
subplot(5,1,3), plot(T,vCB_fix), legend('V_{ONCB}')
subplot(5,1,4), plot(T,w_fix), legend('w')
subplot(5,1,5), plot(T,gbar_MIS_T), legend('g_m')

figure(11)
plot(vCB_fix,w_fix)
xlabel('V_{ONCB}')
ylabel('w')
title(' ONCB Phase Plane')

end


%%%%%%%%%%%%----------------- Function called in run_RelaxOscl.m -----------------------------%%%%%%%%%%%
% Units: mF, mV,ms,S,ohm,mA, cm

function dy = solve_RelaxOscl_2(t,y1)

global epsilon A vz0 vz1 vz2 vz3 vz4 vz5 w0 gamma v0 exp_amp exp_thresh exp_steep tau_w
global delta_tau v_tau v_xi delta_tauM v_tauM v_xiM
global gjONCB ie ie_tonic ie_pulse t_min t_max
global numAII armsection R capacitance            
global gbar_na gbar_M gbar_A G_gap Gpas epas ena ek      
global vhalfh_na vhalfm_na kh_na km_na vhalfm_M km_M vhalfh_A vhalfm_A vhalfw_A kh_A km_A kw_A f
global mtau_na htau_na mtau_M mtau_A 
global tau capacitanceONCB 
global t_switch t_width gbar_M_IS_1 gbar_M_IS_2 SA_hand

y = reshape(y1,numAII,armsection+2,9);
dy = zeros(numAII,armsection+2,9);

% current injection
if (t>=t_min && t<=t_max)
    ie_t=ie_pulse;
else
    ie_t=ie_tonic;
end


gbar_MIS_t=gbar_M_IS_1+(gbar_M_IS_2-gbar_M_IS_1)*.5*(1+tanh((t-t_switch)/t_width));
gbar_M(1,armsection+2)=gbar_MIS_t * SA_hand;
%%%%%%%%%%%% ---------------------------- Differential Equations ------------------------- %%%%%%%%%%%
for i=1:numAII
    
    for j=1:armsection+2
        if (1==2)
            mtau_Mv=mtau_M*(1+delta_tauM/(1+exp(-(y(i,j,8)-v_tauM)/v_xiM)));
        else
            mtau_Mv=mtau_M;
        end
        dy(i,j,1) = (1/mtau_na)*(1/(1+exp(-(y(i,j,7)-vhalfm_na)/km_na))-y(i,j,1));  % m_na (AII)
        dy(i,j,2) = (1/htau_na)*(1/(1+exp((y(i,j,7)-vhalfh_na)/kh_na))-y(i,j,2));  %h_na (AII)
        dy(i,j,3) = (1/mtau_Mv)*(1/(1 + exp(-(y(i,j,7)-vhalfm_M)/km_M))-y(i,j,3));  %m_M (AII)
        dy(i,j,4) = (1/mtau_A)*(1/(1 + exp(-(y(i,j,7)-vhalfm_A)/km_A))-y(i,j,4));  %m_A (AII)
        dy(i,j,5) = (1/(-20/(1+exp(-(y(i,j,7)+35)/6))+25))*(f*(1/(1 + exp((y(i,j,7)-vhalfh_A)/kh_A)))+(1-f)-y(i,j,5));  %h0_A (AII)
        dy(i,j,6) = (1/min(100,((y(i,j,7)+17)^2/4+26)))*(f*(1/(1 + exp((y(i,j,7)-vhalfh_A)/kh_A)))+(1-f)-y(i,j,6));  %h1_A (AII)
        

        if j==1
            isection = 1/R(j)*(y(i,j+1,7)-y(i,j,7));
        elseif j==armsection+2
            isection = 1/R(j-1)*(y(i,j-1,7)-y(i,j,7));
        else
            isection = 1/R(j)*(y(i,j+1,7)-y(i,j,7))+1/R(j-1)*(y(i,j-1,7)-y(i,j,7));
        end
 
        if numAII == 1
            igap = 0;
        else
            if i==1
                igap = G_gap(i,j)*(y(i+1,j,7)-y(i,j,7));
            elseif i==numAII
                igap = G_gap(i,j)*(y(i-1,j,7)-y(i,j,7));
            else
                igap = G_gap(i,j)*(y(i+1,j,7)-y(i,j,7))+G_gap(i,j)*(y(i-1,j,7)-y(i,j,7));
            end
        end

        im = gbar_na(i,j)*y(i,j,1).^3.*y(i,j,2).*(y(i,j,7) - ena)+gbar_M(i,j)*y(i,j,3)*(y(i,j,7) - ek)+1/(1+exp(-(y(i,j,7)-vhalfw_A)/kw_A))*gbar_A(i,j)*y(i,j,4)*y(i,j,5)*(y(i,j,7) - ek) + (1-1/(1+exp(-(y(i,j,7)-vhalfw_A)/kw_A)))*gbar_A(i,j)*y(i,j,4)*y(i,j,6)*(y(i,j,7) - ek)+Gpas(i,j)*(y(i,j,7)-epas);
        
        if (j==1)
            i_fromONCB = gjONCB*(y(i,j,8)-y(i,j,7));
        else
            i_fromONCB=0;
        end
        
%        dy(i,j,7) = (1/capacitance(i,j))*(-im + isection + igap + i_fromONCB + ie); % V_AII
        dy(i,j,7) = (1/capacitance(i,j))*(-im + isection + igap + i_fromONCB + ie_t); % V_AII
        
        if j == 1
            coupling = 1;   %j=1 is the soma
        else
            coupling = 0;
        end
        
        if (1==1)
            dy(i,j,8) =  coupling*((1/capacitanceONCB) * gjONCB*(y(i,j,7)-y(i,j,8))...
                +(1/tau)*(1/epsilon)*(A*(y(i,j,8)-vz0).*(y(i,j,8)-vz1).*(y(i,j,8)-vz2)./(-y(i,j,8)+vz3)./(y(i,j,8)-vz4) - y(i,j,9) - w0+exp_amp/(1+exp(-(y(i,j,8)-exp_thresh)/exp_steep)))); % v_ONCB
        else
            dy(i,j,8) =  coupling*((1/capacitanceONCB) * gjONCB*(y(i,j,7)-y(i,j,8))...
                +(1/tau)*(1/epsilon)*(A*(y(i,j,8)-vz0).*(y(i,j,8)-vz1).*(y(i,j,8)-vz2).*(y(i,j,8)-vz5)./(-y(i,j,8)+vz3) - y(i,j,9) - w0+exp_amp/(1+exp(-(y(i,j,8)-exp_thresh)/exp_steep)))); % v_ONCB
        end
        if (1==1)
            dy(i,j,9) = coupling * 1/(tau*tau_w*(1+delta_tau/(1+exp(-(y(i,j,8)-v_tau)/v_xi))))*(y(i,j,8) - gamma*y(i,j,9) - v0); %w
        else
            dy(i,j,9) = coupling * (1/tau)*(y(i,j,8) - gamma*y(i,j,9) - v0); %w
        end
        
%         gjONCB
        
    end
end

dy = reshape(dy,numAII*(armsection+2)*9,1);

end