Location: BG_Ks @ ef4b8640984c / Clancy_matlab_parameter_fitting / PSO_GHK_fitting_curve.m

Author:
Shelley Fong <s.fong@auckland.ac.nz>
Date:
2022-04-08 14:02:01+12:00
Desc:
Fitting to kinetic hybrid model of LRd + Kernik
Permanent Source URI:
https://models.physiomeproject.org/workspace/82d/rawfile/ef4b8640984c620ac0534355377346919cffdeb6/Clancy_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;

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

cKo = 4.5;
cKi = 141.2;
cNao = 132;
cNai = 9;
Cai = 6e-5; %millimolar
PNaK = 0.01833; % [=] dimless
Cm = 153400e-9; % Unit microF
SA = 11400e-8; % [=] cm2, using SA(um2) = 0.3*vol (in um3) (Bers)

g_Ks = 0.433; % [=] milliS_per_cm2
G_pmr = g_Ks*SA*(1+(0.6/(1+((3.8e-5/Cai)^1.4)))); % [=] mS

% E_K = R*T/F*log(cKo/cKi_st);
% E_K_st = R*T/F*log(cKo_st/cKi_st);
E_K = (R*T/F)*log((cKo+PNaK*cNao)/(cKi+PNaK*cNai)); % [=] V

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

fitRange = 44:67;
% fitRange = [1:length(I_lin)];
error_func = @(G_GHK) square_error(I_lin(fitRange) - calc_IGHK(G_GHK,V(fitRange),cKi,cKo));

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

options_ps = optimoptions('particleswarm','UseParallel',false,'HybridFcn',@fminunc,'SwarmSize',2000, ...
'FunctionTolerance', 1e-16);
if run_optimisation
    [G_GHK,fval,exitflag,output] = particleswarm(error_func,1,lb,ub,options_ps);
    G_GHK2 = G_GHK*0.5;
    save([storage_dir 'Ks_G_GHK.mat'],'G_GHK2');
else
    load([storage_dir 'Ks_G_GHK.mat']);
end

% problem is that sign of current is switched around -80mV, which is
% resting physiological. Want to minimise difference as much as possible.
% hacked to save the more smaller value.
I_GHK = calc_IGHK(G_GHK,V,cKi,cKo);
I_GHK2 = calc_IGHK(G_GHK*0.5,V,cKi,cKo);
I_GHK3 = calc_IGHK(G_GHK*2,V,cKi,cKo);

if run_optimisation
    save([storage_dir 'Ks_G_GHK.mat'],'G_GHK2');
end

h = figure;
plot(1000*V,1e6*I_lin,'k--',1000*V,1e6*I_GHK,'k','LineWidth',2); hold on,
plot(1000*V,1e6*I_GHK2,'r','LineWidth',2); hold on,
plot(1000*V,1e6*I_GHK3,'b','LineWidth',2);
legend('LRd','BG','BG2','BG3','Location','southeast');
ylabel('Current (nA)');
xlabel('Voltage (mV)');
set(gca,'FontSize',16);
grid on;

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