Location: BG_fast_Na @ 6ddc754dd912 / matlab_parameter_fitting / gate_fitting.m

Author:

Shelley Fong <s.fong@auckland.ac.nz>

Date:

2022-04-22 16:55:28+12:00

Desc:

Adding matlab curve and gate fitting files from Pan

Permanent Source URI:

https://models.physiomeproject.org/workspace/826/rawfile/6ddc754dd9127bbf5edbabb61cd6eeb48875a277/matlab_parameter_fitting/gate_fitting.m

% based on markov model of Pan's Na m channel
% work in SECONDS not ms.

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(-100:1:30); % unit mV
% V(91) = []; % remove discontinuity at V=-30 

%  gate : % units: millivolt, milliseconds

% Xf_inf = 1./(1+exp(-(V-1.5)/16.7));
% tau_Xf = 0.001./((7.19e-5*(V+30)./(1-exp(-0.148.*(V+30))))+(0.000131*(V+30)./(exp(0.0687.*(V+30))-1)));
% alpha_Xf = Xf_inf./tau_Xf;
% beta_Xf = (1./tau_Xf) - alpha_Xf;

x_F = [0.0435,5.7897E-7,-14.589712170199999,20086.650237884372,10.202352845280002,23.94529134653];
xF1=x_F(2); xF2=x_F(3); xF5=x_F(4); xF6=x_F(5);
xF3=xF5*xF1;  xF4=1/((1/xF2)+(1/xF6));
xF_const=x_F(6);

% dimensionless variables
alpha_Xf=xF1.*exp((V)./xF2);
beta_Xf=xF3.*exp((V)./xF4);
Xf_inf=alpha_Xf./(alpha_Xf+beta_Xf);
tau_Xf=((1./(alpha_Xf+beta_Xf))+xF_const);

%% Fit bond graph parameters to model
% params: [kf, zf, kr, zr];
error_func_alpha = @(params) square_error(alpha_Xf - calc_alpha(params,V/1000)); % fit to alpha - k.exp(zFV/RT) : do lsqminerror.
error_func_beta = @(params) square_error(beta_Xf - calc_beta(params,V/1000));
error_func_gss = @(params) square_error(Xf_inf - p2gss(params,V));
error_func_tau = @(params) square_error(tau_Xf- p2tau(params,V));

error_func = @(params) error_func_alpha(params) + error_func_beta(params) + error_func_gss(params) + error_func_tau(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',750);

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

V_full = transpose(-120:1:60); % unit mV
alpha_Xf=xF1.*exp((V_full)./xF2);
beta_Xf=xF3.*exp((V_full)./xF4);
Xf_inf=alpha_Xf./(alpha_Xf+beta_Xf);
tau_Xf=((1./(alpha_Xf+beta_Xf))+xF_const);

g_ss_fit = p2gss(params_vec,V_full);
h1 = figure;
plot(V_full,Xf_inf,'k--','LineWidth',4);
hold on;
plot(V_full,g_ss_fit,'k','LineWidth',4);
xlabel('Voltage (mV)');
ylabel('g_{ss}');
legend('linear','Fitted');
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_full);
h2 = figure;
plot(V_full,tau_Xf,'k--','LineWidth',4);
hold on;
plot(V_full,tau_fit,'k','LineWidth',4);
legend('linear','Fitted');
xlabel('Voltage (mV)');
ylabel('\tau_xs1 (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');

alpha_fit = calc_alpha(params_vec,V_full/1000);
beta_fit = calc_beta(params_vec,V_full/1000);
figure,
subplot(1,2,1)
plot(V_full,alpha_Xf,V_full,alpha_fit);
legend('original','fit');
title('alpha');
subplot(1,2,2)
plot(V_full,beta_Xf,V_full,beta_fit);
legend('original','fit');
title('beta');