Location: BG_K_ATP @ 32803a860ae9 / matlab_parameter_fitting / PSO_GHK_fitting_curve.m

Author:
Shelley Fong <s.fong@auckland.ac.nz>
Date:
2022-03-16 14:51:20+13:00
Desc:
Init
Permanent Source URI:
https://models.physiomeproject.org/workspace/83a/rawfile/32803a860ae9213f139af3d11559db3c09f57cee/matlab_parameter_fitting/PSO_GHK_fitting_curve.m
Compare... Download Source 

clear;
% clc;
% close all;

%% Options
run_optimisation = false;

%% Set up directories
current_dir = pwd;
Idx_backslash = find(current_dir == filesep);
main_dir = current_dir; %(1:Idx_backslash(end));
output_dir = [main_dir '\output' filesep];
storage_dir = current_dir; %[main_dir '\storage' filesep];

%% Define constants - coonsistent with VOLTS.
R = 8.314;
T = 310;
F = 96485;

%% Plot I-V curves
% UNITS:
%     I = G_GHK * [K] [=] mA
%     G_GHK = I/[K] [=] Amp.litre/mol
%     G_lin [=] mS [=] mA/V
    
V = (-110:1:65)/1000;

cKo = 4.5;   % mM
cKi = 141.2;
z = 1;
Cm = 153400e-9; % Unit microF

%G_lin = 0.004*Cm; % g_Cab: milliS_per_microF {init: 0.003016}; 
nicholsarea = 5e-5; % dimensionless
i_K_ATP_on = 1; % dimensionless
hATP = 2; kATP = 0.00025; nATP = 0.24;
ATPi = 3; % mM 

g_K_ATP = i_K_ATP_on*0.000193/nicholsarea; % milliS_per_microF
pATP = 1/(1+((ATPi/kATP) ^ hATP)); % dimensionless
G_lin = Cm*g_K_ATP*pATP*((cKo/4) ^nATP); % unit mS
        
E_K = R*T/(z*F)*log(cKo/cKi);

I_lin = G_lin*(V-E_K); % Unit mA    don't consider gating variable: just finding conductance.

Vstart = 1;
Vend =  length(V);
error_func = @(G_GHK) square_error(I_lin(Vstart:Vend) - calc_IGHK(G_GHK,V(Vstart:Vend),cKi,cKo));

A = [];
b = [];
Aeq = [];
beq = [];
lb = [-Inf];
ub = [Inf];

options_ps = optimoptions('particleswarm','UseParallel',false,'HybridFcn',@fminunc,'SwarmSize',1000, ...
    'FunctionTolerance', 1e-9);

if run_optimisation
    [G_GHK,fval,exitflag,output] = particleswarm(error_func,1,lb,ub,options_ps);
    save([storage_dir '\kATP_G_GHK.mat'],'G_GHK');
else
    load([storage_dir '\kATP_G_GHK.mat']);
end

G_GHK2 = 0.00001;

I_GHK = calc_IGHK(G_GHK,V,cKi,cKo);
I_GHK2 = calc_IGHK(G_GHK2,V,cKi,cKo);

h = figure;
% subplot(1,2,1)
plot(1000*V,1e6*I_lin,'k--',1000*V,1e6*I_GHK);
legend('linear','GHK fit','Location','southeast');
ylabel('Current (nA)');
xlabel('Voltage (mV)');
set(gca,'FontSize',16);

% subplot(1,2,2)
% plot(1000*V,1e6*I_lin,'k--',1000*V, 1e6*I_GHK2);
% legend('LRd','BG_test','Location','southeast');
% ylabel('Current (nA)');
% xlabel('Voltage (mV)');
% set(gca,'FontSize',16);

% print_figure(h,output_dir,'tcc_IV_curve');