Location: Han 2019 - Solving a century-old conundrum underlying cardiac force-length relations @ 902a5740b3e6 / ESZ_PlotAllFigures.m

Author:
Kenneth Tran <k.tran@auckland.ac.nz>
Date:
2019-01-08 16:52:34+13:00
Desc:
Uploading Matlab files.
Permanent Source URI:
https://models.cellml.org/workspace/567/rawfile/902a5740b3e6ac6b210dc436d5f25a82ac90070a/ESZ_PlotAllFigures.m
Compare... Download Source 

% ESZ framework code to plot all model simulation figures in manuscript
% This script calls the function 'createFigure.m' to generate the
% background ESZ plot, on top of which the replications are plotted.
% 13 Sept 2018
% Author: Kenneth Tran

%% Run this section first before running the individual sections for plotting
%  the figures in the manuscript

clear all;

% Define functions here
% Quadratic function for isometric work-loop stress
quadS = @(x,a,b,c) a*x.^2 + b*x + c;    

% Exponetial function for the passive stress
expS = @(x,a,b) a*exp(b*x);             

% Solving the quadratic function for the end-systolic length
quadL = @(a,b,c) (-b+sqrt(b*b-4*a.*c))./(2*a);

% Parameters

fs = 12; % Set figure font size

% Parameters for isometric ESSLR
a_isom = 581.3;
b_isom = -852.63; 
c_isom = 314.41;

% Parameters for passive relation
a_pass = 3.289e-7;
b_pass = 17.202;

% Parameters for work-loop ESSLR at Lo 
a_WL1 = 451.9; 
b_WL1 = -489.96;
c_WL1 = 81.709;

% Parameters for work-loop ESSLR at 0.95 Lo 
a_WL2 = 730.79;
b_WL2 = -1066.87;
c_WL2 = 382.73;

% Parameters for work-loop ESSLR at 0.9 Lo 
a_WL3 = 215.28;
b_WL3 = -177.93;
c_WL3 = 4.3608;

% Parameters for normalised work-loop ESSLR
% Used to construct ESSLRs at any Li
a_WL = 15.0202;
b_WL = -19.8935;
c_WL = 5.8743;

%% %%% Figure 3B
fig3B = createFigure(1);

Li = [0.94 0.905 0.87 0.83 0.94 0.905 0.87 0.83];
AL = [0.8 0.77 0.8 0.88 1 1 1 1];

AL_stress = AL.*quadS(Li,a_isom,b_isom,c_isom); % Scales relative afterload to stress (for plotting)
Les = quadL(a_WL,b_WL,c_WL-AL).*Li; % Calculate end-systolic length given Li and AL

% Fit end-systolic points using 2nd order regression
xline = linspace(min(Les)-0.005,max(Les)+0.005);
pline = polyfit(Les,AL_stress,2);
yline = polyval(pline,xline);
plot(xline,yline,'r:','linewidth',2.5)

% Plot isometric and work-loops
for i = 1:length(Les)
    
    EndDiastolicPassive(i) = expS(Li(i),a_pass,b_pass);
    EndSystolicPassive(i) = expS(Les(i),a_pass,b_pass);
    
    switch i
        case {1,2,3,4}
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k','linewidth',1.5)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k','linewidth',1.5)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k','linewidth',1.5)
        otherwise
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k--','linewidth',1)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k--','linewidth',1)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k--','linewidth',1)  
    end
        
    switch i
        case {1,2,3,4}
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','k','markersize',5)
        otherwise
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','w','markersize',5)
    end
  
end

%% %%% Figure 3D
fig3D = createFigure(1);

Li = [0.9 0.89 0.88 0.87 0.899 0.88];
AL = [0.443 0.49497 0.55617 0.6288 0.6 0.4];

AL_stress = AL.*quadS(Li,a_isom,b_isom,c_isom); % Scales relative afterload to stress (for plotting)
Les = quadL(a_WL,b_WL,c_WL-AL).*Li; % Calculate end-systolic length given Li and AL

% Fit end-systolic points using 2nd order regression
xline = linspace(min(Les)-0.005,max(Les)+0.005);
pline = polyfit(Les,AL_stress,2);
yline = polyval(pline,xline);
plot(xline,yline,'r:','linewidth',2.5)

% Plot isometric and work-loops
for i = 1:length(Les)
    
    EndDiastolicPassive(i) = expS(Li(i),a_pass,b_pass);
    EndSystolicPassive(i) = expS(Les(i),a_pass,b_pass);
    
    switch i
        case {1,2,3,4}
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k','linewidth',1)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k','linewidth',1)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k','linewidth',1)
        otherwise
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k:','linewidth',1.5)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k:','linewidth',1.5)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k:','linewidth',1.5)
    end
        
    switch i
        case {1,2,3,4}
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','k','markersize',5)
        otherwise
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','w','markersize',5)
    end
  
end

%% %%% Figure 3F
fig3F = createFigure(1);

Li = [0.9 0.885 0.847 0.86 0.878];
AL = [0.7 0.65 1 1 1];

AL_stress = AL.*quadS(Li,a_isom,b_isom,c_isom); % Scales relative afterload to stress (for plotting)
Les = quadL(a_WL,b_WL,c_WL-AL).*Li; % Calculate end-systolic length given Li and AL

% Fit end-systolic points using 2nd order regression
xline = linspace(min(Les)-0.005,max(Les)+0.005);
pline = polyfit(Les,AL_stress,2);
yline = polyval(pline,xline);
plot(xline,yline,'r:','linewidth',2.5)

% Plot isometric and work-loops
for i = 1:length(Les)
    
    EndDiastolicPassive(i) = expS(Li(i),a_pass,b_pass);
    EndSystolicPassive(i) = expS(Les(i),a_pass,b_pass);
    
    switch i
        case {1,2}
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k','linewidth',1)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k','linewidth',1)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k','linewidth',1)
        otherwise
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k--','linewidth',1)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k--','linewidth',1)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k--','linewidth',1)
    end
        
    switch i
        case {1,2}
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','k','markersize',5)
        otherwise
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','w','markersize',5)
    end
  
end

%% %%% Figure 3H
fig3H = createFigure(1);

Li = [0.9 0.9 0.9 0.9 0.901 0.8808 0.864 0.8475 0.8285];
AL = [0.8 0.63 0.47 0.3 1 1 1 1 1];

AL_stress = AL.*quadS(Li,a_isom,b_isom,c_isom); % Scales relative afterload to stress (for plotting)
Les = quadL(a_WL,b_WL,c_WL-AL).*Li; % Calculate end-systolic length given Li and AL

% Fit end-systolic points using 2nd order regression
xline = linspace(min(Les)-0.005,max(Les)+0.005);
pline = polyfit(Les,AL_stress,2);
yline = polyval(pline,xline);
plot(xline,yline,'r:','linewidth',2.5)

% Plot isometric and work-loops
for i = 1:length(Les)
    
    EndDiastolicPassive(i) = expS(Li(i),a_pass,b_pass);
    EndSystolicPassive(i) = expS(Les(i),a_pass,b_pass);
    
    switch i
        case {1,2,3,4}
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k','linewidth',1)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k','linewidth',1)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k','linewidth',1)
        otherwise
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k--','linewidth',1)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k--','linewidth',1)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k--','linewidth',1)
    end
        
    switch i
        case {1,2,3,4}
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','k','markersize',5)
        otherwise
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','w','markersize',5)
    end
  
end

%% %%% Figure 4C
fig4C = createFigure(0);

Li = [0.93 0.92 0.895 0.875 0.93 0.928 0.927 0.925 0.92 0.91 0.90 0.895 0.875];
AL = [0.76 0.68 0.68 0.58 1 1 1 1 1 1 1 1 1];

AL_stress = AL.*quadS(Li,a_isom,b_isom,c_isom); % Scales relative afterload to stress (for plotting)
Les = quadL(a_WL,b_WL,c_WL-AL).*Li; % Calculate end-systolic length given Li and AL

% Plot isometric and work-loops
for i = 1:length(Les)
    
    EndDiastolicPassive(i) = expS(Li(i),a_pass,b_pass);
    EndSystolicPassive(i) = expS(Les(i),a_pass,b_pass);
    
    switch i
        case {1,2,3,4}
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k','linewidth',1.5)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k','linewidth',1.5)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k','linewidth',1.5)
        case {5,9,12,13}
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k--','linewidth',0.5)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k--','linewidth',0.5)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k--','linewidth',0.5)               
    end
        
    switch i
        case {1,2,3,4}
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','k','markersize',5)
        otherwise
            plot(Les(i),AL_stress(i),'kx','markersize',6)
    end
  
end

%% %%% Figure 4D
fig4D = createFigure(0);

Li = [0.945 0.92 0.88 0.85 0.945 0.928 0.913 0.901 0.92 0.87 0.86 0.88 0.85];
AL = [0.4 0.45 0.45 0.6 1 1 1 1 1 1 1 1 1];

AL_stress = AL.*quadS(Li,a_isom,b_isom,c_isom); % Scales relative afterload to stress (for plotting)
Les = quadL(a_WL,b_WL,c_WL-AL).*Li; % Calculate end-systolic length given Li and AL

% Plot isometric and work-loops
for i = 1:length(Les)
    
    EndDiastolicPassive(i) = expS(Li(i),a_pass,b_pass);
    EndSystolicPassive(i) = expS(Les(i),a_pass,b_pass);
    
    switch i
        case {1,2,3,4}
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k','linewidth',1.5)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k','linewidth',1.5)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k','linewidth',1.5)
        case {5,9,12,13}
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k--','linewidth',0.5)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k--','linewidth',0.5)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k--','linewidth',0.5)               
    end
        
    switch i
        case {1,2,3,4}
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','k','markersize',5)
        otherwise
            plot(Les(i),AL_stress(i),'kx','markersize',6)
    end
  
end

%% %%% Figure 5B
fig5B = createFigure(0,[0.8,0.85,0.9],[0,10,20],[0.76,0.92],[0,20]);

Li = [0.90 0.90 0.90 0.90 0.90 0.90 0.89 0.89 0.89 0.89 0.89 0.89 0.8805 0.8805 0.8805 0.8805 0.8805];
AL = [1 0.8 0.6 0.43 0.25 0.1 1 0.8939 0.6706 0.4805 0.2794 0.1117 1 0.7491 0.5369 0.3122 0.1248];

AL_stress = AL.*quadS(Li,a_isom,b_isom,c_isom); % Scales relative afterload to stress (for plotting)
Les = quadL(a_WL,b_WL,c_WL-AL).*Li; % Calculate end-systolic length given Li and AL

% Fit end-systolic points using 2nd order regression
xline = linspace(min(Les),max(Les));
pline = polyfit(Les,AL_stress,2);
yline = polyval(pline,xline);
plot(xline,yline,'r','linewidth',1.5)

% Plot isometric and work-loops
for i = 1:length(Les)
    
    EndDiastolicPassive(i) = expS(Li(i),a_pass,b_pass);
    EndSystolicPassive(i) = expS(Les(i),a_pass,b_pass);
    
    switch i
        case {1,7,13}
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k--','linewidth',1)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k--','linewidth',1)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k--','linewidth',1)
        otherwise
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k','linewidth',1)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k','linewidth',1)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k','linewidth',1)
    end
        
    switch i
        case {1,2,3,4,5,6}
            plot(Les(i),AL_stress(i),'ks','markerfacecolor','k','markersize',6)
        case {7,8,9,10,11,12}
            plot(Les(i),AL_stress(i),'ks','markerfacecolor','w','markersize',6)
        case {13,14,15,16,17}
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','w','markersize',6)
    end
    
    switch i
        case {1}
            plot(Les(i),AL_stress(i),'ks','markerfacecolor','k','markersize',6)
        case {7}
            plot(Les(i),AL_stress(i),'ks','markerfacecolor','w','markersize',6)
        case {13}
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','w','markersize',6)
    end

end

%% %%% Figure 5C
fig5C = createFigure(0,[0.8,0.85,0.9],[0,10,20],[0.76,0.92],[0,20]);

Li = [0.90 0.90 0.90 0.90 0.90 0.90 0.89 0.89 0.89 0.89 0.89 0.89 0.8805 0.8805 0.8805 0.8805 0.8805];
AL = [1 0.8 0.6 0.43 0.25 0.1 1 0.8939 0.6706 0.4805 0.2794 0.1117 1 0.7491 0.5369 0.3122 0.1248];

AL_stress = AL.*quadS(Li,a_isom,b_isom,c_isom); % Scales relative afterload to stress (for plotting)
Les = quadL(a_WL,b_WL,c_WL-AL).*Li; % Calculate end-systolic length given Li and AL

% Fit end-systolic points using 2nd order regression
xline = linspace(min(Les)+0.015,max(Les));
pline = polyfit([Les(1) Les(2) Les(3) Les(4) Les(5) Les(6)],[AL_stress(1) AL_stress(2) AL_stress(3) AL_stress(4) AL_stress(5) AL_stress(6)],2);
yline = polyval(pline,xline);
plot(xline,yline,'r:','linewidth',2)

xline = linspace(min(Les)+0.004,max(Les));
pline = polyfit([Les(7) Les(8) Les(9) Les(10) Les(11) Les(12)],[AL_stress(7) AL_stress(8) AL_stress(9) AL_stress(10) AL_stress(11) AL_stress(12)],2);
yline = polyval(pline,xline);
plot(xline,yline,'r:','linewidth',2)

xline = linspace(min(Les)-0.005,max(Les)-0.02);
pline = polyfit([Les(13) Les(14) Les(15) Les(16) Les(17)],[AL_stress(13) AL_stress(14) AL_stress(15) AL_stress(16) AL_stress(17)],2);
yline = polyval(pline,xline);
plot(xline,yline,'r:','linewidth',2)


% Plot isometric and work-loops
for i = 1:length(Les)
    
    EndDiastolicPassive(i) = expS(Li(i),a_pass,b_pass);
    EndSystolicPassive(i) = expS(Les(i),a_pass,b_pass);
    
    switch i
        case {1,7,13}
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k--','linewidth',1)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k--','linewidth',1)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k--','linewidth',1)
        otherwise
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k','linewidth',1)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k','linewidth',1)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k','linewidth',1)
    end
        
    switch i
        case {1,2,3,4,5,6}
            plot(Les(i),AL_stress(i),'ks','markerfacecolor','k','markersize',6)
        case {7,8,9,10,11,12}
            plot(Les(i),AL_stress(i),'ks','markerfacecolor','w','markersize',6)
        case {13,14,15,16,17}
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','w','markersize',6)
    end
    
    switch i
        case {1}
            plot(Les(i),AL_stress(i),'ks','markerfacecolor','k','markersize',6)
        case {7}
            plot(Les(i),AL_stress(i),'ks','markerfacecolor','w','markersize',6)
        case {13}
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','w','markersize',6)
    end

end

%% %%% Figure 6B
fig6B = createFigure(0);

Li = 0.93.*[1 1 1 1 1 1 0.97];
AL = [1 0.87 0.73 0.53 0.32 0.12 1];

AL_stress = AL.*quadS(Li,a_isom,b_isom,c_isom); % Scales relative afterload to stress (for plotting)
Les = quadL(a_WL,b_WL,c_WL-AL).*Li; % Calculate end-systolic length given Li and AL

% Fit end-systolic points using 2nd order regression
xline = linspace(min(Les),max(Les));
pline = polyfit(Les(1:6),AL_stress(1:6),2);
yline = polyval(pline,xline);
plot(xline,yline,'r:','linewidth',2.5)

% Plot isometric and work-loops
for i = 1:length(Les)
    
    EndDiastolicPassive(i) = expS(Li(i),a_pass,b_pass);
    EndSystolicPassive(i) = expS(Les(i),a_pass,b_pass);
    
    switch i
        case {1,7}
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k','linewidth',1)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k','linewidth',1)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k','linewidth',1)
        otherwise
            plot([Li(i) Li(i)],[EndDiastolicPassive(i) AL_stress(i)],'k','linewidth',1)
            plot([Li(i) Les(i)],[AL_stress(i) AL_stress(i)],'k','linewidth',1)
            plot([Les(i) Les(i)],[EndSystolicPassive(i) AL_stress(i)],'k','linewidth',1)
    end
        
    switch i
        case {1,7}
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','w','markersize',5)
        otherwise
            plot(Les(i),AL_stress(i),'ko','markerfacecolor','k','markersize',5)
    end
  
end