% Run this script for simulating the response of  the 1:1 BK-CaV complex to
% an AP generated by a neuron model. In particular, simulate the responses
% obtained from the 6-state ODE model (Eqs. S6-S11), 
% the simplified Hodgkin-Huxley-type model (Eq. 23 coupled with Eq. 20)
% and the further simplified model assuming instantaneous activation (Eq. 29
% with k=1, m_{CaV}=m_{CaV,\infty}, coupled with Eq. 20),
% as shown in Figure S5. Also, compare the results with 
% those obtained from the 70-state Markov chain model (Cox, 2014).

clear all 
clc
close all

plot_cox=1; % if one plot the results obtained from the Cox model

% define global parameters  
global k1_0 K1  k2_0 K2 n1 n2 alpha0 beta0 V_0  ca_c  n_c  k_minus k_plus N_tot beta1 alpha1 beta1_0 alpha1_0 rho
% BK parameters
bk_par=[  1.1093   3.3206  2.3298    0.0223    1.6012   16.5793    0.1000    0.4614];

k1_0=bk_par(1);

k2_0=bk_par(2);

K1=bk_par(6);

K2=bk_par(7);

n1=bk_par(3);

n2=bk_par(8);

beta0=bk_par(4);

alpha0=-bk_par(5)*beta0;

ca_c=0.2; % background Ca2+ concentration
V_0=-80; %initial voltage value

n_c=1; % number of calcium channels (CaVs) coupled with BK
% ca_o= -6/20*V+19;

%CaV and BK channnels
n_cav_tot=1000;
n_bk_tot=1000;
% constants for Ampere to nA and pA
conv_na=1e-12*1e-3*1e9;
conv_pa=1e-12*1e-3*1e12;

N_tot=1; % the probabilities sum is 1 (N_tot=cx+cy+ox+oy+bx+by)

% CaV parameters
cv_par=[1.2979    1.0665   -0.0639    0.0703    0.3090]; 

alpha1_0=cv_par(1);

beta1_0=cv_par(2);

alpha1=cv_par(3);

beta1=cv_par(4);

rho=cv_par(5);

% load the data for the action potential
load('data_AP');

        
% window time for AP
tv=time_AP;
time_in=tv(1);
time_f=tv(end);

% action potential 
V_in=v_AP;

%%% compute calcium level at 7nm of the pore (cav sensor) and at 13nm (bk sensor)
D_ca=250; %microm^2 s^-1
F=9.6485*10^4; %C mol^-1
conv_F=10^(-15);  % mol a M/microm^3
conv_microM=10^6; % M to microM
r_bk=13*10^(-3);% microm
r_ca=7*10^(-3); %microm
k_B=500; %microM^-1 * s^-1
B_tot=30; %microM
Eca=60; %mV
E_K=-96;
conv_V=10^(-3); %mV to V
g_ca=2.8; %pS
g_bk=100;%pS
conv_S=10^(-12); %pS to S
i_ca=abs((V_in-Eca))*conv_V*g_ca*conv_S; % C/sec

ca_o_rca=i_ca/(8*pi*D_ca*F*conv_F*r_ca)*exp(-r_ca/(sqrt(D_ca/k_B/B_tot)))*conv_microM;
ca_o=i_ca/(8*pi*D_ca*F*conv_F*r_bk)*exp(-r_bk/(sqrt(D_ca/k_B/B_tot)))*conv_microM;
% ca_o= -6/20*V+19;

% CaV parameters
alpha=alpha1_0.*exp(-alpha1.*V_in); % alpha in the paper
beta=beta1_0.*exp(-beta1.*V_in);
beta_scale=(beta+alpha)* rho; % beta in the paper
m_inf=alpha./(alpha+beta_scale); % m_{CaV,\infty}

k1=k1_0.*exp(-alpha0.*V_in); % w^{+} defined in the paper
k2=k2_0.*exp(-beta0.*V_in); % w^{-} defined in the paper

fco_Ca_c=(ca_c^n1)/(ca_c^n1+K1^n1); % f^{+}(Ca_c) defined in the paper
fco_Ca_o=(ca_o.^n1)./(ca_o.^n1+K1^n1); % f^{+}(Ca_o) defined in the paper

foc_Ca_c=(K2^n2)/(K2^n2+ca_c^n2); % f^{-}(Ca_c) defined in the paper
foc_Ca_o=(K2^n2)./(K2^n2+ca_o.^n2); % f^{-}(Ca_o) defined in the paper


k1_ca_o=k1.*fco_Ca_o; % k^{+}_{o} defined in the paper
k1_ca_c=k1.*fco_Ca_c; % k^{+}_{c} defined in the paper
k2_ca_o=k2.*foc_Ca_o; % k^{-}_{o} defined in the paper
k2_ca_c=k2.*foc_Ca_c; % k^{-}_{c} defined in the paper

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BK-CaV model with n_c=1 non-inactivated CaV%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tspan=time_in:1e-2:time_f;
x_0=zeros(1, 2*(n_c+1)) ;
x_0(1)=1;

% 4-state ODE model for the BK activation
[T, Y]=ode45('bk_n_cav_states_model_scale',tspan, x_0,[],tv,V_in);
bk_y=Y(:,n_c+2:end);
y_sum=sum(bk_y,2); % BK open prob, p_y, with non-inactivated CaV


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BK-CaV  model approximation with n_c=1 non-inactivated CaV%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yk1_0=0;
yk2_0=[0,1];

% 1-ODE model for the BK activation (Eq.29 with k=1) with m_CaV= m_{CaV,\infty}
[T1, y1qssminf]=ode45('bk_cav_1state_model_scale_qss_minf',tspan, yk1_0,[],tv,V_in);

% 1-ODE model for the BK activation (Eq. 23) and 1-ODE model for the CaV activation
tspan2=time_in:1e-2:time_f;

[T1a, Y2]=ode45('bk_cav_2state_model_scale_qss',tspan2, yk2_0,[],tv,V_in);

yk2=Y2(:,1); % m_{BK}
c2=Y2(:,2);  % c in the paper for CaV


%%%%%%%%%%%%%%%%%%%%%%
%%%%%% CaV model %%%%%%%
%%%%%%%%%%%%%%%%%%%%%%
tspan1=time_in:1e-2:time_f;
x_10=[1 0];
k_minus=2*10^(-3); % gamma in the paper
k_plus=2.5*10^(-3); % delta/Ca_{CaV} in the paper

% 2-state ODE model for CaV
[T2c, Y2c]=ode45('cav_2states_model_scale',tspan1, x_10,[],tv,V_in);
c=Y2c(:,1);
o=Y2c(:,2);
b=1-c-o;
h_ca=1-b;

% 1-state ODE model for CaV
[T1, b1]=ode45('cav_1state_h_scale',tspan1, 0,[],tv,V_in);

h1=1-b1; % h for CaV with 1-state ODE model

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% BK-CaV model with n_c=1 inactivating CaV %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tspan6= time_in:1e-2:time_f;
x_0=[0, 0,0,0,0 ];
% 6-state ODE model for 1:1 BK-CaV
[T6, Y6]=ode45('bk_cav_6states_model_scale',tspan6, x_0,[],tv,V_in);

oy=Y6(:,1);
cy=Y6(:,2);
ox=Y6(:,3);
bx=Y6(:,4);
by=Y6(:,5);
cx=N_tot-oy-cy-ox-bx-by;

y=oy+cy+by; % p_y, BK open prob with the 6-state ODE model
y_a_tot_1=y1qssminf.*h1; % BK open prob with 1-state ODE for m_{BK} (Eq. 29 with k=1), m_{CaV}=m_{CaV,\infty}, and 1-state ODE model for h
y_a_tot_2=yk2.*h1; % BK open prob with 1-state ODE for m_{BK} (Eq. 23), 1-state ODE for m_{CaV}, and 1-state ODE model for h

% for figure size
wl=200;
wh=200;

% CaV open prob with instantaneous activation
minf_h_CaV= m_inf.*h1;

t_sh=1;
x_inf=1-t_sh;
x_sup=6-t_sh;

f=figure(1);
hold on
grid on
plot(T-t_sh,o,'-b','linewidth',1);
plot(T-t_sh,minf_h_CaV,'-.g','linewidth',1);
xlim([0 40])
set(f,'Position',[10 10 wl wh])
xlabel('time [ms]');
ylabel('CaV open prob.')


f=figure(2);
hold on
grid on
plot(T-t_sh,conv_pa*n_cav_tot*g_ca*o.*(V_in-Eca),'-b','linewidth',1);
plot(T-t_sh,conv_pa*n_cav_tot*g_ca*minf_h_CaV.*(V_in-Eca),'-.g','linewidth',1);
xlim([x_inf x_sup])
set(f,'Position',[10 10 wl wh])
xlabel('time [ms]');
ylabel('I_{CaV} [pA]')
title('Figure S5B')


f=figure (3);
hold on
grid on
plot(T6-t_sh,y,'-b','linewidth',1); %
plot(T-t_sh,y_a_tot_2,'--r','linewidth',1);
plot(T-t_sh,y_a_tot_1,'-.g','linewidth',1);
set(f,'Position',[10 10 wl wh])
xlim([x_inf x_sup])
xlabel('time [ms]');
ylabel('BK open prob.')


f=figure (4);
hold on
grid on
title('Figure S5C')
plot(T6-t_sh,conv_na*n_bk_tot*g_bk*y.*(V_in-E_K),'-b','linewidth',1); %
plot(T-t_sh,conv_na*n_bk_tot*g_bk*y_a_tot_2.*(V_in-E_K),'--r','linewidth',1);
plot(T-t_sh,conv_na*n_bk_tot*g_bk*y_a_tot_1.*(V_in-E_K),'-.g','linewidth',1);
set(f,'Position',[10 10 wl wh])
xlim([x_inf x_sup])
xlabel('time [ms]');
ylabel('I_{BK} [nA]')


f=figure(5);
title('Figure S5A')
hold on
grid on
plot(T-t_sh,V_in,'-b','linewidth',1);
xlim([x_inf x_sup])
ylim([-80 41])
set(f,'Position',[10 10 wl wh])
xlabel('time [ms]');
ylabel('AP [mV]')

if plot_cox
     
    load ('MC_sim/Cox_model_AP.mat')
    figure(1);
    hold on
    grid on
    plot(t-t_sh,open_cav_cox_AP,':','color',[0.5,0.5,0.5],'linewidth',1)
    xlim([x_inf x_sup])

    figure(3);
    hold on
    grid on
    plot(t-t_sh,open_bk_cox_AP,':','color',[0.5,0.5,0.5],'linewidth',1)
    
    xlim([x_inf x_sup])

    figure(2);
    hold on
    grid on
    plot(t-t_sh,conv_pa*n_cav_tot*g_ca*open_cav_cox_AP.*(V_in-Eca)',':','color',[0.5,0.5,0.5],'linewidth',1)
    xlim([x_inf x_sup])

    figure(4);
    hold on
    grid on
    plot(t-t_sh,conv_na*n_bk_tot*g_bk*open_bk_cox_AP.*(V_in-E_K)',':','color',[0.5,0.5,0.5],'linewidth',1)
    xlim([x_inf x_sup])
       
end