Location: BG_funny @ b2e59dcb052c / Kernik_matlab_parameter_fitting / PSO_GHK_fitting_curve.m

Author:
Shelley Fong <sfon036@UoA.auckland.ac.nz>
Date:
2022-07-06 16:24:56+12:00
Desc:
Updating to LRd 34.4 pL cell parameters
Permanent Source URI:
https://models.cellml.org/workspace/838/rawfile/b2e59dcb052c70b8a3675561fd41a23c48846604/Kernik_matlab_parameter_fitting/PSO_GHK_fitting_curve.m
Compare... Download Source 

clear;
% clc;
% close all;

%% Options
run_optimisation = true;

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

%% Define constants
R = 8.314;
T = 310;
F = 96485;

%% choose to plot either Na or K portions of the current
Na = true;

%% Plot I-V curves
% UNITS:
%     I = G_GHK * K [=] mA
%     G_GHK = I/K [=] Amp.litre/mol
%     G_PMR [=] mS [=] mA/V
    
V = (-120:1:60)/1000;
fitRange = [1:length(V)];

cKo = 4.5;
cKi = 141.2;
cNao = 132;
cNai = 9;
Cai = 6e-5; %millimolar
PNaK = 0.01833; % [=] dimless
Cm = 153400e-9; % Unit microF

E_Na = R*T/F*log(cNao/cNai);
E_K = R*T/F*log(cKo/cKi);

G_0 = 0.0435; % nS_per_pF = mS/uF
NatoK_ratio = 0.491; %Verkerk et al. 2013
Na_frac = NatoK_ratio./(NatoK_ratio+1);
if Na
    G_f = Na_frac*G_0*Cm;
    I_lin = G_f*(V-E_Na); % Unit mA    don't consider gating variable: just finding conductance.
    error_func = @(G_GHK) square_error(I_lin(fitRange) - calc_IGHK(G_GHK,V(fitRange),cNai,cNao));
else
    G_f = (1-Na_frac)*G_0*Cm;
    I_lin = G_f*(V-E_K); % Unit mA    don't consider gating variable: just finding conductance.
    error_func = @(G_GHK) square_error(I_lin(fitRange) - calc_IGHK(G_GHK,V(fitRange),cKi,cKo));
end


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

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

if run_optimisation
    [G_GHK,fval,exitflag,output] = particleswarm(error_func,1,lb,ub,options_ps);
    if Na
        save([storage_dir 'fNa_G_GHK.mat'],'G_GHK');
    else
        save([storage_dir 'fK_G_GHK.mat'],'G_GHK');
    end
else
    if Na
        load([storage_dir 'fNa_G_GHK.mat']);
    else
        load([storage_dir 'fK_G_GHK.mat']);
    end
end
if Na
    I_GHK = calc_IGHK(G_GHK,V,cNai,cNao);
else
    I_GHK = calc_IGHK(G_GHK,V,cKi,cKo);
end

h = figure;
plot(1000*V,1e6*I_lin,'k--',1000*V,1e6*I_GHK,'k','LineWidth',2);
legend('LRd','BG','Location','southeast');
ylabel('Current (nA)');
xlabel('Voltage (mV)');
if Na
    title('Na')
else
    title('K')
end
set(gca,'FontSize',16);
grid on