- Author:
- Shelley Fong <sfon036@UoA.auckland.ac.nz>
- Date:
- 2025-08-04 14:38:27+12:00
- Desc:
- uptdating diag with gating
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/d54/rawfile/c75f331b160c1ce02438ef644c9fb3c90fa7d02e/parameter_finder/find_G_GHK.m
% curve fitting for i_ns_Ca, which consists of 3 currents (each carrying
% Ca, Na, K)
% 3 GHK fits must be done to find 3 GHK conductances
clear;
% clc;
% close all;
%% Options
run_optimisation = true;
%% 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 = [main_dir '\storage' filesep];
%% Define constants - consistent with VOLTS.
R = 8.314;
T = 310;
F = 96485;
%% choose which ion to plot current of
ion = 'K'; %Ca, Na, K
%% 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;
V(find(V == 0,1)) = [];
V = V';
cCao = 2.5;
cCai = 0.000116; %millimolar
cNao = 130;
cNai = 4;
cKo = 6;
cKi = 140;
% Tong
A = 7.07e-12*1e4; % cm2
Cm_per = 1; % uF_per_cm2
Cm = Cm_per*A*1e6; % pF
g_kinetic_per = 0.0123; % nS_per_pF
g_kinetic = g_kinetic_per*Cm;
G = g_kinetic*1e-6; % Unit mA/V = mS
% LRd
% Cm = 153400e-9; % Unit microF
% SA = 11400e-6; % [=] mm2, using SA(um2) = 0.3*vol (in um3) (Bers)
% P_ns_Ca = 1.75e-7; % [=] cm_per_second
% K_m_ns_Ca = 0.0012; % [=] mM
% gamma = 0.75; % [=] dimensionless
if false
% offset, the same as LRd94 code
EnsCa = R*T/F*log((cKo+cNao)/(cKi+cNai));
Vns = V - EnsCa;
else
Vns = V;
end
%
% % the following currents in units [=] uA_per_mm2
% I_ns_Na = P_ns_Ca*Vns*(F^2)/(R*T).*(gamma*cNai*exp(1*Vns*F/(R*T))-gamma*cNao)./(exp(Vns*F/(R*T))-1);
% I_ns_K = P_ns_Ca*Vns*(F^2)/(R*T).*(gamma*cKi*exp(1*Vns*F/(R*T))-gamma*cKo)./(exp(Vns*F/(R*T))-1);
% i_ns_Na = I_ns_Na./(1+(K_m_ns_Ca/cCai ^ 3));
% i_ns_K = I_ns_K./(1+(K_m_ns_Ca/cCai ^ 3));
% i_ns_Ca = i_ns_Na+i_ns_K;
V = (-120:1:60)/1000;
E_Ca = R*T/F*log(cCao/cCai);
E_Na = R*T/F*log(cNao/cNai);
E_K = R*T/F*log(cKo/cKi);
vFRT = V*F/(R*T);
PnsK = 1.3;
PnsNa = 0.9;
PnsCa = 0.89;
PnsCs = 1.0;
E_ns = (R*T/F)*log((PnsK*cKo+PnsNa*cNao+4*PnsCa*cCao./(1.0+exp(vFRT)))./(PnsK*cKi+PnsNa*cNai+4*PnsCa*cCai./(1.0+exp(vFRT))));
% i_ns_Na = G*(V-E_Na); % Unit mA
% i_ns_K = G*(V-E_K); % Unit mA
% i_ns_Ca = G*(V-E_Ca); % Unit mA
i_ns_Na = G*(V-E_ns); % Unit mA
i_ns_K = G*(V-E_ns); % Unit mA
i_ns_Ca = G*(V-E_ns); % Unit mA
Vstart = 1;
Vend = round(length(V)*1);
switch ion
case 'Ca'
z = 2;
error_func = @(G_GHK) square_error(i_ns_Ca(Vstart:Vend) - calc_IGHK(G_GHK,V(Vstart:Vend),cCai,cCao, z));
mat_name = '\nsCa_G_GHK.mat';
case 'Na'
z = 1;
error_func = @(G_GHK) square_error(i_ns_Na(Vstart:Vend) - calc_IGHK(G_GHK,V(Vstart:Vend),cNai,cNao, z));
mat_name = '\nsNa_G_GHK.mat';
case 'K'
z = 1;
error_func = @(G_GHK) square_error(i_ns_K(Vstart:Vend) - calc_IGHK(G_GHK,V(Vstart:Vend),cKi,cKo, z));
mat_name = '\nsK_G_GHK.mat';
otherwise
error('not an ion')
end
A = [];
b = [];
Aeq = [];
beq = [];
lb = [-Inf];
ub = [Inf];
options_ps = optimoptions('particleswarm','UseParallel',false,'HybridFcn',@fminunc,'SwarmSize',1500, ...
'FunctionTolerance', 1e-16);
if run_optimisation
[G_GHK,fval,exitflag,output] = particleswarm(error_func,1,lb,ub,options_ps);
save([storage_dir mat_name],'G_GHK');
else
load([storage_dir mat_name]);
end
disp(G_GHK);
h = figure;
% subplot(1,2,1)
if strcmp(ion,'Ca')
I_GHK = calc_IGHK(G_GHK,V,cCai,cCao,z);
% I_GHK2 = calc_IGHK(G_GHK*0.5,V,cCai,cCao,z);
% I_GHK3 = calc_IGHK(G_GHK*2,V,cCai,cCao,z);
plot(1000*V,1e6*i_ns_Ca,'k--',1000*V,1e6*I_GHK);
% hold on;
% plot(1000*V,1e6*I_GHK2,'r');
% plot(1000*V,1e6*I_GHK3,'b');
elseif strcmp(ion,'Na')
I_GHK = calc_IGHK(G_GHK,V,cNai,cNao,z);
plot(1000*V,1e6*i_ns_Na,'k--',1000*V,1e6*I_GHK);
else
I_GHK = calc_IGHK(G_GHK,V,cKi,cKo,z);
plot(1000*V,1e6*i_ns_K,'k--',1000*V,1e6*I_GHK);
end
legend('LRd','GHK fit','Location','southeast');
ylabel('Current (nA)');
xlabel('Voltage (mV)');
title(ion)
set(gca,'FontSize',16);
grid on
