Location: BG_TCC @ 06f090addb49 / Kernik_matlab_parameter_fitting / PSO_GHK_fitting_curve.m

Author:
Shelley Fong <sfon036@UoA.auckland.ac.nz>
Date:
2024-11-13 10:26:46+13:00
Desc:
Updating cellml and sedml
Permanent Source URI:
https://models.cellml.org/workspace/831/rawfile/06f090addb4968d12531197d463345e0b2012d2a/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_bactcclash = find(current_dir == filesep);
main_dir = current_dir; %(1:Idx_bactcclash(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 = (-150:1:60)/1000;

cCao = 1.8; %0.2454;
cCai = 0.0003; % avg of LRd conc %0.00012; %6e-5;
Cm = 153400e-9; % Unit microF
z = 2;

G_pmr = 0.185*Cm; %Cm*0.433*(1+0.6/(1+((3.8e-5/Cai)^1.4))); 


% test with Kr values
if false
    cCao = 5.4;
    cCai = 104.75;
    z = 1;
%     Cm = 60e-6; % Kernik  % 153400e-9; % Unit microF
    G_pmr = 0.218025*sqrt(cCao/5.4)*Cm; 
end

E_Ca = R*T/(z*F)*log(cCao/cCai);

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

Vstart = 1; %31;
Vend = length(V); %101;
error_func = @(G_GHK) square_error(I_lin(Vstart:Vend) - calc_IGHK(G_GHK,V(Vstart:Vend),cCai,cCao,z));

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

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

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

G_GHK2 = 0.00001;

I_GHK = calc_IGHK(G_GHK,V,cCai,cCao,z);
I_GHK2 = calc_IGHK(G_GHK2,V,cCai,cCao,z);

h = figure;
% subplot(1,2,1)
plot(1000*V,1e6*I_lin,'k--',1000*V,1e6*I_GHK);
legend('LRd','BG','Location','southeast');
ylabel('Current (nA)');
xlabel('Voltage (mV)');
set(gca,'FontSize',16);
grid on
% 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');