Location: BG_TCC @ bc202e51bde6 / Kernik_matlab_parameter_fitting / d_gate_fitting.m

Author:

Shelley Fong <sfon036@UoA.auckland.ac.nz>

Date:

2022-07-06 16:26:28+12:00

Desc:

Updating to LRd 34.4 pL cell parameters

Permanent Source URI:

https://models.fieldml.org/workspace/831/rawfile/bc202e51bde6be8ddf88d35d908ef1465dcf067d/Kernik_matlab_parameter_fitting/d_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(-120:1:80); % unit mV

%  gate : % units: millivolt, ms (same as Pan)
d_inf = 1.0./(1.0+exp((V+26.3)/-6.0));
tau_d = 1.0./(1.068*exp((V+26.3)/30.0)+1.068*exp((V+26.3)/-30.0));

alpha_d = d_inf./tau_d;
beta_d = (1./tau_d) - alpha_d;

%% Fit bond graph parameters to model
% params: [kf, zf, kr, zr];
% Kf corresponds to k_zdv20: 1*alpha_d
error_func_alpha = @(params) square_error(alpha_d - calc_alpha(params,V/1000)); % fit to alpha - k.exp(zFV/RT) : do lsqminerror.
error_func_deta = @(params) square_error(beta_d - calc_beta(params,V/1000));
error_func_gss = @(params) square_error(d_inf - p2gss(params,V));
error_func_tau = @(params) square_error(tau_d- p2tau(params,V));

error_func = @(params) error_func_alpha(params) + error_func_deta(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,'FunctionTolerance',1e-14);

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

%% PLOT
% plot over a physiological V range
if false
    V2 = transpose(-120:1:50);
else
    V2 = V;
end

g_ss_fit = p2gss(params_vec,V2);
h1 = figure;
plot(V2,d_inf,'k.');
hold on;
plot(V2,g_ss_fit,'k');
xlabel('Voltage (mV)');
ylabel('m_{ss}');
legend('original','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,V2);
h2 = figure;
plot(V2,tau_d,'k.');
hold on;
plot(V2,tau_fit,'k');
legend('original','Fitted');
xlabel('Voltage (mV)');
ylabel('\tau_d (ms)');
set(gca,'FontSize',28);
xlim([-120 60]);
set(gca,'XTick',-120:30:60);
xticklabels({-120,'',-60,'',0,'',60});
set(gca,'LineWidth',3);
set(gca,'xgrid','on');

alpha_fit = calc_alpha(params_vec,V2/1000);
beta_fit = calc_beta(params_vec,V2/1000);
figure,
subplot(1,2,1)
plot(V2,alpha_d,V2,alpha_fit);
legend('original','fit');
subplot(1,2,2)
plot(V2,beta_d,V2,beta_fit);
legend('original','fit');

% print_figure(h1,output_dir,'d_inf');
% print_figure(h2,output_dir,'tau_d');