- Author:
- Shelley Fong <sfon036@UoA.auckland.ac.nz>
- Date:
- 2024-11-12 11:16:41+13:00
- Desc:
- Updating cellml and sedml
- Permanent Source URI:
- https://models.cellml.org/workspace/82d/rawfile/6ca9cbf971a30a2b5d74e93b0f3cdc816b065680/Kernik_matlab_parameter_fitting/xs_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 = false;
%% 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:60); % unit mV
% V(91) = []; % remove discontinuity at V=-30
% gate : % units: millivolt, milliseconds
% xs_infinity = 1./(1+exp(-(V-1.5)/16.7));
% tau_xs = 0.001./((7.19e-5*(V+30)./(1-exp(-0.148.*(V+30))))+(0.000131*(V+30)./(exp(0.0687.*(V+30))-1)));
% alpha_xs = xs_infinity./tau_xs;
% beta_xs = (1./tau_xs) - alpha_xs;
alpha_xs = 0.0011656*exp(V/66726.84);
beta_xs = 0.0003269*exp(V/-18.872);
tau_offset = 0;
tau_xs = 1./(alpha_xs + beta_xs) + tau_offset;
xs_infinity = alpha_xs./(alpha_xs + beta_xs); % Pan [=] dimensionless
%% Fit bond graph parameters to model
% params: [kf, zf, kr, zr];
% Kf corresponds to k_zdv20: 1*alpha_xs
error_func_alpha = @(params) square_error(alpha_xs - calc_alpha(params,V/1000)); % fit to alpha - k.exp(zFV/RT) : do lsqminerror.
error_func_beta = @(params) square_error(beta_xs - calc_beta(params,V/1000));
error_func_gss = @(params) square_error(xs_infinity - p2gss(params,V));
error_func_tau = @(params) square_error(tau_xs- 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',2000);
options_ps = optimoptions('particleswarm','UseParallel',false,'HybridFcn',@fmincon,'SwarmSize',1550,...
'FunctionTolerance',1e-15);
if run_optimisation
[params_vec,fval,exitflag,output] = particleswarm(error_func,4,lb,ub,options_ps);
save([storage_dir 'ks_xs_parameters.mat'],'params_vec');
else
load([storage_dir 'ks_xs_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,xs_infinity,'kx','LineWidth',4);
hold on;
plot(V,g_ss_fit,'k','LineWidth',4);
xlabel('Voltage (mV)');
ylabel('m_{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_xs,'kx','LineWidth',4);
hold on;
plot(V,tau_fit,'k','LineWidth',4);
% legend('Luo and Rudy','Fitted');
xlabel('Voltage (mV)');
ylabel('\tau_xs (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/1000);
beta_fit = calc_beta(params_vec,V/1000);
figure,
subplot(1,2,1)
plot(V,alpha_xs,'x',V,alpha_fit);
legend('original','fit');
title('alpha');
subplot(1,2,2)
plot(V,beta_xs,'x',V,beta_fit);
legend('original','fit');
title('beta');
% print_figure(h1,output_dir,'xs_infinity');
% print_figure(h2,output_dir,'tau_xs');
