- Author:
- Kenneth Tran <k.tran@auckland.ac.nz>
- Date:
- 2017-08-31 14:58:36+12:00
- Desc:
- Adding MatLab code for users to reproduce simulation figures in the J Physiol publication: DOI: 10.1113/JP274680.
- Permanent Source URI:
- https://models.cellml.org/workspace/4a2/rawfile/eb49874a619c53bd0b587ea6854313751e76edd3/Main_noPassive.m
% Simulating mechano-energetics of isometric and work-loops contractions
% This code reproduces the NO PASSIVE FORCE simulations in the JPhysiol paper:
% "Experimental and modelling evidence of shortening heat in cardiac muscle"
% Kenneth Tran
% August 2017
% The simulations in this script take about 50s to solve
% on a laptop (i7 6650U CPU @2.20 GHZ)
%%
clear;
% SL for isometric contractions
preload_SL_iso = [1.635 1.935 2.035 2.11 2.175 2.23 2.3];
Afterload = [1]; % Not used in isometric simulations - a dummy value
kxb = 89.69; % Scaling constant to scale normalised force to units of kPa
G_ATP = -58; % Free energy of ATP hydrolysis (kJ/mol)
passive = 0; % Passive present = 1; Passive absent = 0;
% Set the twitch protocol
protocol = 'Isometric';
% Creating matrices for storing variables from ISOMETRIC contractions
p = 800;
Force_iso_vec = nan*ones(p, length(preload_SL_iso));
T_iso_vec = nan*ones(p, length(preload_SL_iso));
ATP_iso_vec = nan*ones(p, length(preload_SL_iso));
SL_iso_vec = nan*ones(p, length(preload_SL_iso));
for i=1:length(preload_SL_iso)
loop = 0; % No-loop simulation
T_loop = [2000 2000]; % Not used in non-loop simulations - these are dummy values
[T_iso Y_iso dy_iso F_total_iso ATPase_iso] = XBSolve(loop,T_loop,preload_SL_iso(i),Afterload(1),protocol,passive);
SL_iso = Y_iso(:,9);
% Storing the vectors
T_iso_vec(1:length(T_iso),i) = T_iso;
Force_iso_vec(1:length(F_total_iso),i) = kxb*F_total_iso;
SL_iso_vec(1:length(SL_iso),i) = SL_iso;
ATP_iso_vec(1:length(ATPase_iso),i) = ATPase_iso;
% Finding the index for t<800ms to avoid numerical spike at end
j = 0;
find = 1;
while(find)
j = j+1;
if (T_iso(j)>800)
find = 0;
end
end
% Integrating to compute total ATP consumed in a twitch
ATPase_iso_per_beat(i) = trapz(T_iso(1:j), ATPase_iso(1:j));
% Isometric twitch Peak Force
% Maximum Peak Force (i.e. PF for SL = 2.3) is used to normalise afterloads and stresses
PF(i) = max(F_total_iso);
end
% Total enthalpy
enthalpy_IS = -ATPase_iso_per_beat*G_ATP;
% Max enthalpy - used to normalise energy outputs
ME = max(enthalpy_IS);
disp('Isometric simulations completed.');
%% Workloop contraction.
% Each set of workloops at a given preload takes approx 20 sec
% on a laptop (i7 6650U CPU @2.20 GHZ)
% String for naming the initial preloads for the workloop contractions
SL_string = {'NoPassive_SL23'};
% Peak isometric force for each preload
% Used to scale the work-loop afterloads
SL_PF = [0.6690];
% Initial preloads for the workloop contractions
preload_SL = [2.3];
% Normalised afterloads for work-loop contractions
% Note its multiplication by SL_PF when it is used in the code
Afterload = [0.02 0.05 0.1 0.15 0.2 0.35 0.5 0.65 0.8 0.99];
kxb = 89.69; % Scaling constant to scale normalised force to units of kPa
G_ATP = -58; % Free energy of ATP hydrolysis (kJ/mol)
% Creating vectors for storing variables for WORKLOOP contractions
p = 800;
T_vec = nan*ones(p, length(Afterload));
Force_vec = nan*ones(p, length(Afterload));
SL_vec = nan*ones(p, length(Afterload));
ATP_vec = nan*ones(p, length(Afterload));
% Set the contraction protocol
protocol = 'Workloop';
for k=1:length(SL_string)
for i=1:length(Afterload)
loop = 0; % No-loop simulation
T_loop = [2000 2000]; % Not used in non-loop simulations - these are dummy values
[T_noL Y_noL dy_noL F_total_noL ATPase_noL] = XBSolve(loop,T_loop,preload_SL(k),SL_PF(k)*Afterload(i),protocol,passive);
SL_noL = Y_noL(:,9);
% Find the minimum SL (SL_min) from the above simulation
% This is used to hold the SL fixed at SL_min until T = 900
[Y I] = min(SL_noL);
T_start = T_noL(I);
T_stop = 900;
T_loop = [T_start T_stop]; %Used to set start time for loop
loop = 1; % Simulation with loop
[T_L Y_L dy_L F_total_L ATPase_L] = XBSolve(loop,T_loop,preload_SL(k),SL_PF(k)*Afterload(i),protocol,passive);
SL_L = Y_L(:,9);
% Restretch phase
SL_WL(i) = SL_L(end); % SL at end of workloop
n = 20;
T_stretch = linspace(901,1000,n);
SL_stretch = linspace(SL_WL(i),preload_SL(k),n);
for j=1:length(SL_stretch)
Force_stretch(j) = passiveForces(SL_stretch(j),passive);
end
% Storing the vectors
T_vec(1:length(T_L),i) = T_L;
T_vec(length(T_L)+1:length(T_L)+20,i)= T_stretch;
Force_vec(1:length(F_total_L),i) = kxb*F_total_L;
Force_vec(length(F_total_L)+1:length(F_total_L)+20,i) = kxb*Force_stretch;
SL_vec(1:length(SL_L),i) = SL_L;
SL_vec(length(SL_L)+1:length(SL_L)+20,i) = SL_stretch;
ATP_vec(1:length(ATPase_L),i) = ATPase_L;
% Finding the index for t<800ms to avoid numerical spike at end
j = 0;
find = 1;
while(find)
j = j+1;
if (T_L(j)>800)
find = 0;
end
end
% Integrating to calculate total ATP consumption and work done
ATPase_per_beat(i) = trapz(T_L(1:j), ATPase_L(1:j));
Work(i) = trapz([SL_L;SL_stretch']/max(preload_SL(k)),kxb*[F_total_L Force_stretch]);
end
enthalpy_WL = -ATPase_per_beat*G_ATP;
work = -Work; % Work is negative for contraction
heat_WL = enthalpy_WL - work;
% Interpolate the isometric enthalpy-stress curve at equivalent
% afterloads in the work-loop curve
enthalpy_IS_interp = interp1(PF/max(PF),enthalpy_IS,SL_PF(k)*Afterload/max(PF));
SH = heat_WL - enthalpy_IS_interp;
SH(end) = 0;
% Outputting the data
eval([SL_string{k},'_T_vec = T_vec;']);
eval([SL_string{k},'_SL_vec = SL_vec;']);
eval([SL_string{k},'_Force_vec = Force_vec;']);
eval([SL_string{k},'_AF = SL_PF(k)*Afterload;']);
eval([SL_string{k},'_enthalpy_IS = enthalpy_IS;']);
eval([SL_string{k},'_enthalpy_WL = enthalpy_WL;']);
eval([SL_string{k},'_work = work;']);
eval([SL_string{k},'_ATP_vec = ATP_vec;']);
eval([SL_string{k},'_heat_WL = heat_WL;']);
eval([SL_string{k},'_enthalpy_IS_interp = enthalpy_IS_interp;']);
eval([SL_string{k},'_SH = SH;']);
disp([SL_string{k},' work-loops completed.']);
end
% Run the PlotFigures_noPassive.m script to
% generate figures from this simulation
