Location: The cardiac Na+ /K+ ATPase: An updated, thermodynamically consistent model @ 6baaa1d0bfa1 / Code / MATLAB / code / CellML_model_comparison.m

Author:
Michael Pan <mpan92@gmail.com>
Date:
2020-08-27 08:39:13+10:00
Desc:
Fixed units
Permanent Source URI:
https://models.physiomeproject.org/workspace/4ce/rawfile/6baaa1d0bfa1a0772c316ab59631b619df2601f4/Code/MATLAB/code/CellML_model_comparison.m
Compare... Download Source 

% This script loads the results from running the kinetic and bond graph
% models in CellML, and plots a comparison (Figure 5)

clear;
clc;
close all;

%% Directory names
current_dir = cd;
Idx_backslash = find(current_dir == filesep);
main_dir = current_dir(1:Idx_backslash(end));
data_dir = [main_dir 'data' filesep];
output_dir = [main_dir 'output' filesep];
storage_dir = [main_dir 'storage' filesep];
  
%% Plot comparison of kinetic and bond graph models
load([storage_dir 'kinetic_fitting_results.mat']);
% Constants
struct_input.T = 310;

% Species concentrations
struct_input.Ki = 0;
struct_input.Ke = 5.4;
struct_input.Nai = 50;
struct_input.Nae = 150;
struct_input.MgATP = 10;
struct_input.MgADP = 0;
struct_input.Pi_total = 0;
struct_input.H_conc = 10^(-4.4);
struct_input.V = -0.08;

Vm_vec = (-120:1:60)/1000;

h = figure;
hold on;
box on;

% Nae = 150;
% cycling_velocity_vec = zeros(length(Vm_vec),1);
% for i_Vm = 1:length(Vm_vec)
%     V_m = Vm_vec(i_Vm);
%     struct_input.Nae = Nae;
%     struct_input.V = V_m;
% 
%     cycling_velocity_vec(i_Vm) = NaK_fitting_vcyc(params_vec,struct_input);
% end
% plot(1000*Vm_vec,cycling_velocity_vec,'LineWidth',2);

CellML_kinetic_data = csvread([storage_dir 'Terkildsen_NaK_kinetic_Data.csv'],1,0);
V_kinetic_data = CellML_kinetic_data(:,3);
vcyc_kinetic_data = CellML_kinetic_data(:,2);
plot(V_kinetic_data,vcyc_kinetic_data,'k','LineWidth',3);

CellML_BG_data = csvread([storage_dir 'Terkildsen_NaK_BG_Data.csv'],1,0);
V_BG_data = 1000*CellML_BG_data(:,2);
vcyc_BG_data = CellML_BG_data(:,3);
plot(V_BG_data,vcyc_BG_data,'r--','LineWidth',3);

ylim([10 60]);
xlabel('Membrane voltage (mV)');
ylabel('Cycling velocity (s^{-1})');
lgd = legend('Kinetic','BG','Location','southeast');
set(gca,'FontSize',18);

print_figure(h,output_dir,'CellML_comp');