Location: BG_TCC @ b34e6db3f3f0 / Clancy_matlab_parameter_fitting / g_gate_fitting.m

Author:
Shelley Fong <s.fong@auckland.ac.nz>
Date:
2022-04-08 16:52:55+12:00
Desc:
Reformulate gate to be same as Kernik
Permanent Source URI:
https://models.cellml.org/workspace/831/rawfile/b34e6db3f3f05814a462fa69110445cf287fbbb1/Clancy_matlab_parameter_fitting/g_gate_fitting.m
Compare... Download Source 

% based on Na_h gate by Pan.

clear;
% clc;
% close all;

%% Options
run_optimisation = true;

%% Set directories
current_dir = pwd;
main_dir = current_dir; %
data_dir = [main_dir '\data' filesep];
output_dir = [main_dir '\output' filesep];
storage_dir = [main_dir '\storage' filesep];

%% Steady-state gating parameters and time constants
V = transpose(-120:4:60);

g_inf = 1./(1+exp((V+60)/5.6));
tau_g = max(-0.000875*V+0.012, 0.012);

alpha_g = g_inf./tau_g;
beta_g = (1./tau_g) - alpha_g;

%% Fit bond graph parameters to model
% params: [kf, zf, kr, zr];
error_func_alpha = @(params) square_error(alpha_g - calc_alpha(params,V/1000));
error_func_beta = @(params) square_error(beta_g - calc_beta(params,V/1000));
error_func_gss = @(params) square_error(g_inf - p2gss(params,V));
error_func_tau = @(params) square_error(tau_g- p2tau(params,V));

% error_func = @(params) error_func_alpha(params) + error_func_beta(params) + error_func_gss(params) + error_func_tau(params);
% error_func = @(params) error_func_gss(params) + error_func_tau(params);
error_func = @(params) error_func_alpha(params) + error_func_beta(params);


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

options_unc = optimoptions('fminunc','MaxFunEvals',10000);
options_ps = optimoptions('particleswarm','UseParallel',false,'HybridFcn',@fmincon,'SwarmSize',1750);

if run_optimisation
    [params_vec,fval,exitflag,output] = particleswarm(error_func,4,lb,ub,options_ps);
    save([storage_dir 'TCC_g_parameters.mat'],'params_vec');
else
    load([storage_dir 'TCC_g_parameters.mat']);
end
% [params_vec,fval,exitflag,output,grad,hessian] = fminunc(error_func,params_vec,options_unc);

g_ss_fit = p2gss(params_vec,V);
h1 = figure;
plot(V,g_inf,'k.','LineWidth',4);
hold on;
plot(V,g_ss_fit,'k','LineWidth',4);
legend('original','Fitted');
xlabel('Voltage (mV)');
ylabel('h_{ss}');
set(gca,'FontSize',28);
xlim([-120 60]);
set(gca,'XTick',-120:30:60);
set(gca,'YTick',0:0.2:1);
xticklabels({-120,'',-60,'',0,'',60});
yticklabels({0,'','','','',1});
set(gca,'LineWidth',3);
grid on;

tau_fit = p2tau(params_vec,V);
h2 = figure;
plot(V,tau_g,'k.','LineWidth',4);
hold on;
plot(V,tau_fit,'k','LineWidth',4);
legend('original','Fitted');
xlabel('Voltage (mV)');
ylabel('\tau_g (ms)');
set(gca,'FontSize',28);
xlim([-120 60]);
set(gca,'XTick',-120:30:60);
% set(gca,'YTick',0:0.2:1);
xticklabels({-120,'',-60,'',0,'',60});
% yticklabels({0,'','','','',1});
set(gca,'LineWidth',3);
set(gca,'xgrid','on');

print_figure(h1,output_dir,'g_inf');
print_figure(h2,output_dir,'tau_g');

alpha_fit = calc_alpha(params_vec,V/1000);
beta_fit = calc_beta(params_vec,V/1000);
figure,
subplot(1,2,1)
plot(V,alpha_g,V,alpha_fit);
legend('original','fit');
subplot(1,2,2)
plot(V,beta_g,V,beta_fit);
legend('original','fit');