- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-09-03 12:15:33+12:00
- Desc:
- Init
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/6fa/rawfile/e459ccc5beeea2965725c7eea6409c97c3bbd563/matlab/main.m
% coupled signaling/EC for adult rat ventricular myocytes
%
% FIGURE REPRODUCTION MODE
% Run this script and follow prompts regarding which figure in the primary paper to reproduce.
% Reference:
% Jeffrey J. Saucerman, Laurence L. Brunton, Anushka p. Michailova, and Andrew D. McCulloch
% "Modeling beta-adrenergic Control of cardiac myocyte contractility in silico", J. Biol. Chem., Vol 278: 47977-48003
%
% Copyright (2004) The Regents of the University of California
% All Rights Reserved
%
% Last modified: 3 September 2021
% Implemented by: Shelley Fong
%
clear;
close all;
tic
% -----------------------------------
% checkFigure 2a 2b 2c 2d 2e 2f 2gh 3a 3b 4a 4b 5 6a 6b
checkFigure = input('Enter figure number in format ia: ', 's');
% -----------------------------------
if strcmp(checkFigure ,'2g') || strcmp(checkFigure ,'2h')
checkFigure = '2gh';
end
fprintf('Fig %s\n',checkFigure)
saucerman_jbc2003(checkFigure);
toc
function saucerman_jbc2003(checkFigure)
extracted_figure_path = 'figures/paperExtractedCSV/';
isoproterenol = 0.01;
protocol = 0;
cAMP_holding = repmat([0],50,1);
PKI = repmat([0],50,1);
expr_mult = repmat([0],50,1);
CTX = repmat([0],50,1);
vPLBKO = repmat([0],50,1);
IBMX = repmat([0],50,1);
data_fracPLBp = [];
data_fig4 = [];
data_cAMP_normal = [];
data_cAMP_ACmult = [];
data_cAMP_B1ARmult = [];
data_cAMP_KGsadiv = [];
data_cAMP_Control = [];
data_cAMP_PKA_KO = [];
data_cAMP_PKI_KO = [];
data_CaTransient = [];
data_CaTransient_PLBKO = [];
% --------------- EXPERIMENTAL PROTOCOLS ---------------
N = 1e3;
switch checkFigure
case '2a'
tspan = [0,18*60];
numloops0 = 1;
numloops = numloops0;
isoproterenol = 0.01;
case '2b'
tspan = linspace(0,5*60,N);
numloops0 = 8; % different isoproterenol. *2 for IBMX.
numloops = numloops0*2;
isoproterenol = logspace(-3,0.45,numloops0); % solution diverges at isoproterenol > 10^0.1 uM
data_cAMP_ss_IBMX0 = zeros(numloops0,1);
data_cAMP_ss_IBMX1 = zeros(numloops0,1);
IBMX0 = [0, 100]; % uM % 270
% repeat loops to run with IBMX
IBMX = [IBMX0(1)*ones(numloops0,1); IBMX0(2)*ones(numloops0,1)];
isoproterenol = [isoproterenol, isoproterenol];
case {'2c','2d'}
tspan = linspace(0,60,N);
numloops0 = 6; % num trials with different [cAMP] *2 for PKI.
isoproterenol = ones(numloops0*2,1)*isoproterenol;
cAMP_holding = logspace(-1.9, 2, numloops0);
PKI0 = [0e-15, 0.052];
PKI = [PKI0(1)*ones(numloops0,1); PKI0(2)*ones(numloops0,1)];
lenPKI = length(PKI0);
numloops = numloops0*lenPKI;
cAMP_holding = [cAMP_holding, cAMP_holding];
data_cAMP_ss0 = zeros(numloops0*2,2); % []
data_PKA0 = zeros(numloops0*2,1); % []
data_cAMP_ss1 = zeros(numloops0*2,2); % []
data_PKA1 = zeros(numloops0*2,1); % []
case '2e'
tspan = linspace(0, 60*30,N); % 30 minutes
numloops0 = 20;
isoproterenol = logspace(-6,0,numloops0);
numloops = numloops0;
case '2f'
tspan = [0:0.001:60];
numloops0 = 2; % without/with isoproterenol
numloops = numloops0;
isoproterenol = [0, 1];
protocol = 2;
case '2gh'
tspan = [0:0.001:60];
numloops0 = 2; % without/with isoproterenol
numloops = numloops0;
isoproterenol = [0, 1];
protocol = 1;
case {'3a', '3b'}
numloops0 = 8; % vary isoproterenol conc
if strcmp(checkFigure, '3a')
tspan = linspace(0,7*60,N);
isoproterenol = logspace(-3,0,numloops0);
expr_mult0 = [1,3.0, 3.0, 3.0]; % [Control, ACmult, B1ARmult, cAMP_mult/K_Gsa] [1,3,3,3]
IBMX0 = 1e2;
else
tspan = linspace(0,6*60,N);
isoproterenol = logspace(-4,-0.5,numloops0);
expr_mult0 = [1,1.5, 1.8, 1.9]; % [Control, ACmult, B1ARmult, cAMP_mult/K_Gsa]
IBMX0 = 0;
end
IBMX = repmat(IBMX0,numloops0*4,1);
lenx = length(expr_mult0);
numloops = lenx*numloops0;
expr_mult = [expr_mult0(1)*ones(numloops0,1); expr_mult0(2)*ones(numloops0,1); expr_mult0(3)*ones(numloops0,1); expr_mult0(4)*ones(numloops0,1)];
isoproterenol = repmat(isoproterenol, 1, lenx);
case '4a'
tspan = linspace(0, 3*60);
numloops0 = 3;
numloops = numloops0;
isoproterenol = [1,0,1];
IBMX0 = 2.5e2;
IBMX = repmat(IBMX0,numloops0,1);
CTX = [0,1,1];
case '4b'
tspan = linspace(0, 3*60);
numloops0 = 2;
isoproterenol = [1e-4, 1];
CTX = [0,0,1,1];
IBMX0 = 2.5e2;
numloops = numloops0*2;
isoproterenol = repmat(isoproterenol,1,2);
IBMX = repmat(IBMX0,numloops,1);
data_Gsa_gtp_AC = zeros(numloops,1);
case '5'
tspan = linspace(0, 20*60,N); % 20*60
numtrials = 3;
numloops0 = 8; % different isoproterenol
isoproterenol = logspace(-4,0,numloops0);
isoproterenol = repmat(isoproterenol, 1, numtrials);
% Control PKA_KO PKI_KO
numloops = numloops0*numtrials;
case '6a'
tspan = [0:0.001:60];
numloops0 = 3; % normal, LCCphos only, PLB phos only.
protocol = 1;
numloops = numloops0;
isoproterenol = ones(numloops, 1);
case '6b'
tspan = [0:0.001:60];
numloops0 = 4; % normal, Iso, PLBKO, Iso+PLBKO
protocol = 1;
assert(numloops0 == 4, 'numloops must be 4')
numloops = numloops0;
isoproterenol = [0,1,0,1];
vPLBKO = [0,0,1,1];
otherwise
error('no figure specified');
end
if false && (~strcmp(checkFigure,'2f')&&(~strcmp(checkFigure,'2gh') && (~strcmp(checkFigure,'6a') && ~strcmp(checkFigure,'6b'))))
for it = 1:length(tspan)
tspan(it) = tspan(it)*0.001;
end
end
Cm = 324e-6; % uF: 207 or 324 pF
% ----- Signaling model parameters -------
% b-AR/Gs module
p_orig(1) = 0; %isoproterenol(n); % Ltotmax [uM]
p_orig(2) = 0.0132; % sumb1AR [uM]
p_orig(3) = 3.83; % Gstot [uM]
p_orig(4) = 0.285; % Kl [uM]
p_orig(5) = 0.062; % Kr [uM]
p_orig(6) = 33.0; % Kc [uM]
p_orig(7) = 1.1e-3; % k_barkp [1/sec]
p_orig(8) = 2.2e-3; % k_barkm [1/sec]
p_orig(9) = 3.6e-3; % k_pkap [1/sec/uM]
p_orig(10) = 2.2e-3; % k_pkam [1/sec]
p_orig(11) = 16.0; % k_gact [1/sec]
p_orig(12) = 0.8; % k_hyd [1/sec]
p_orig(13) = 1.21e3; % k_reassoc [1/sec/uM]
% cAMP module
p_orig(14) = 49.7e-3;% AC_tot [uM]
p_orig(15) = 5.0e3; % ATP [uM]
p_orig(16) = 38.9e-3;% PDEtot [uM]
p_orig(17) = 0.0; % IBMXtot [uM]
p_orig(18) = 0.0; % Fsktot [uM] (10 uM when used)
p_orig(19) = 0.2; % k_ac_basal[1/sec]
p_orig(20) = 8.5; % k_ac_gsa [1/sec]
p_orig(21) = 7.3; % k_ac_fsk [1/sec]
p_orig(22) = 5.0; % k_pde [1/sec]
p_orig(23) = 1.03e3; % Km_basal [uM]
p_orig(24) = 315.0; % Km_gsa [uM]
p_orig(25) = 860.0; % Km_fsk [uM]
p_orig(26) = 1.3; % Km_pde [uM]
p_orig(27) = 0.4; % Kgsa [uM]
p_orig(28) = 44.0; % Kfsk [uM]
p_orig(29) = 30.0; % Ki_ibmx [uM]
% PKA module
p_orig(30) = 0.59; % PKAItot [uM]
p_orig(31) = 0.025; % PKAIItot [uM]
p_orig(32) = 0.18; % PKItot [uM] 0.18
p_orig(33) = 9.14; % Ka [uM]
p_orig(34) = 1.64; % Kb [uM]
p_orig(35) = 4.375; % Kd [uM]
p_orig(36) = 0.2e-3; % Ki_pki [uM]
% PLB module
p_orig(37) = 10; % epsilon [none]
p_orig(38) = 106; % PLBtot [uM]
p_orig(39) = 0.89; % PP1tot [uM]
p_orig(40) = 0.3; % Inhib1tot [uM]
p_orig(41) = 54; % k_pka_plb [1/sec]
p_orig(42) = 21; % Km_pka_plb [uM]
p_orig(43) = 8.5; % k_pp1_plb [1/sec]
p_orig(44) = 7.0; % Km_pp1_plb [uM]
p_orig(45) = 60; % k_pka_i1 [1/sec]
p_orig(46) = 1.0; % Km_pka_i1 [uM]
p_orig(47) = 14.0; % Vmax_pp2a_i1 [uM/sec]
p_orig(48) = 1.0; % Km_pp2a_i1 [uM]
p_orig(49) = 1.0e-3; % Ki_inhib1 [uM]
% LCC module
p_orig(50) = 25e-3; % LCCtot [uM]
p_orig(51) = 25e-3; % PP1lcctot [uM]
p_orig(52) = 25e-3; % PP2Alcctot [uM]
p_orig(53) = 54; % k_pka_lcc [1/sec]
p_orig(54) = 21; % Km_pka_lcc [uM]
p_orig(55) = 8.52; % k_pp1_lcc [1/sec]
p_orig(56) = 3; % Km_pp1_lcc [uM]
p_orig(57) = 10.1; % k_pp2a_lcc [1/sec]
p_orig(58) = 3; % Km_pp2a_lcc [uM]
% ---- EC Coupling model parameters ------
% universal parameters
p_orig(59) = 20.8e-6; % Vmyo [uL]
p_orig(60) = 9.88e-7; % Vnsr [uL]
p_orig(61) = 9.3e-8; % Vjsr [uL]
p_orig(62) = 1.534e-4; % ACap [cm^2] with C = 1 uF/cm^2
p_orig(63) = 310; % Temp [K]
% extracellular concentrations
p_orig(64) = 140; % Extracellular Na [mM]
p_orig(65) = 5.4; % Extracellular K [mM]
p_orig(66) = 1.8; % Extracellular Ca [mM]
% current conductances
p_orig(67) = 8.0; % G_Na [mS/uF]
p_orig(68) = 0.35; % G_to [mS/uF]
p_orig(69) = 0.07; % G_ss [mS/uF]
p_orig(70) = 0.24; % G_kibar [mS/uF]
p_orig(71) = 0.008; % G_kp [mS/uF]
% I_Ca parameters
p_orig(72) = 300; % f [1/sec]
p_orig(73) = 2e3; % g [1/sec]
p_orig(74) = 5187.5; % gammao [1/sec/mM]
p_orig(75) = 10; % omega [1/sec]
p_orig(76) = 5.823e-9*3.0; % pCa [cm/sec]
p_orig(77) = 1.078e-11*3.0; % pK [cm/sec]
p_orig(78) = 3e5; % Nlcc [#/cell]
p_orig(79) = -0.458; % I_Ca05 [uA/uF]
% pumps and background currents
p_orig(80) = 1483; % k_NaCa [uA/uF]
p_orig(81) = 87.5; % Km_Na [mM]
p_orig(82) = 1.38; % Km_Ca [mM]
p_orig(83) = 0.1; % k_sat [none]
p_orig(84) = 0.35; % eta [none]
p_orig(85) = 1.1; % ibarnak [uA/uF]
p_orig(86) = 10; % Km_Nai [mM]
p_orig(87) = 1.5; % Km_Ko [mM]
p_orig(88) = 1.15; % ibarpca [uA/uF]
p_orig(89) = 0.5e-3; % Km_pca [mM]
p_orig(90) = 2.8e-3; % G_Cab [uA/uF]
p_orig(91) = 1.18e-3; % G_Nab [uA/uF]
p_orig(92) = 0; % Pns
p_orig(93) = 1.2e-3; % Km_ns [mM]
% Calcium handling parameters
p_orig(94) = 4.7; % I_upbar [mM/sec]
p_orig(95) = 3e-4; % Km_up [mM]
p_orig(96) = 15; % nsrbar [mM]
p_orig(97) = 2e-3; % tauon [sec]
p_orig(98) = 2e-3; % tauoff [sec]
p_orig(99) = 60e3; % gmaxrel [mM/sec]
p_orig(100) = 0.18e-3;% dcaith [mM]
p_orig(101) = 0.8e-3; % Km_rel [mM]
p_orig(102) = 8.75; % CSQNth [mM]
p_orig(103) = 15; % CSQNbar [mM]
p_orig(104) = 0.8; % Km_csqn [mM]
p_orig(105) = 5.7e-4; % tau_tr [sec]
p_orig(106) = 0.07; % TRPNbar [mM]
p_orig(107) = 0.05; % CMDNbar [mM]
p_orig(108) = 0.07; % INDObar [mM]
p_orig(109) = 0.5128e-3; % Km_trpn [mM]
p_orig(110) = 2.38e-3; % Km_cmdn [mM]
p_orig(111) = 8.44e-4; % Km_indo [mM]
% Initial conditions and mass matrix
% b-AR/Gsa module
% 1 2 3 4 5 6 7 8 9
% L R Gs b1ARtot b1ARd b1ARp Gsagtptot Gsagdp Gsbg
y10=[0.010; 0.01079;3.829; 0.01205;0.0; 1.154e-3;0.02505; 6.446e-4;0.02569];
m1 =[0, 0, 0, 1, 1, 1, 1, 1, 1];
% cAMP module and PKA module
% 10 11 12 13 14 15
% Gsa_gtp Fsk AC PDE IBMX cAMPtot
y20=[0.02241;0; 0.04706;0.0389; 0; 0.8453];
m2 =[0, 0, 0, 0, 0, 1];
% PKA module
% 16 17 18
% cAMP PKACI PKACII
y30=[0.2268;0.05868;8.278e-3];
m3 =[0, 0, 0];
% PLB
% 19 20 21 22
% PLBp Inhib1ptot Inhib1p PP1
y40=[4.105; 0.0526; 6.3e-5; 0.838];
m4 =[1, 1, 0, 0];
% LCC module
% 23 24
% LCCap LCCbp
y50=[5.103e-3; 5.841e-3];
m5 =[1, 1];
% Gating variables
% 25 26 27 28 29 30 31 32 33
% m h j v w x y z rto sto ssto rss sss
y60=[1.4e-3;0.99; 0.99; 0.0; 0.0; 0.13; 0.96; 0.92; 1.4e-3; 1.0; 0.613; 198e-3; 0.43];
m6 =[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1];
% Intracellular concentrations/ Membrane voltage
% 38 39 40 41 42 43 44
% Ca_nsr Ca_jsr Nai Ki Cai Vm trel
y70=[1.92; 1.92; 16; 145; 1.58e-4;-85.66; 0.9];
m7 =[1, 1, 1, 1, 1, 1, 1];
% Put everything together
y0_orig = [y10;y20;y30;y40;y50;y60;y70];
M = diag([m1,m2,m3,m4,m5,m6,m7]);
% folder for timecourse plots
temp_path = 'figures/temp/';
if ~exist(temp_path,'dir')
mkdir(temp_path)
else
% clear folder
dinfo = dir(temp_path);
dinfo([dinfo.isdir]) = [];
filenames = fullfile(temp_path, {dinfo.name});
if length(filenames)
delete( filenames{:} )
end
end
% convert cAMP concentration (uM) to pmol/mg
% volume of one myocyte = 36.8 pL (Bers)
Vmyo = 36.8e-12;
fudgeFactor = 6.4; %6*15/14; % to fit with data
cellsPermg = 25635.76702; %Mendez and Keys 1961 (from myocyte density)
cAMP_mult = Vmyo*cellsPermg*1e6*fudgeFactor;
% fprintf('fudge factor = %d for [cAMPtot]\n',fudgeFactor);
% disp(numloops)
parfor (n=1:numloops)
% for n=1:numloops
p = p_orig;
y0 = y0_orig;
fprintf('Loop %d/%d \n', n,numloops);
CTXn = 0;
plot_cAMP = false;
fig6p = '';
fracPLBo = 0.9613;
PLBKO = 0;
switch checkFigure
case {'2a', '2b', '2e'}
p(1) = isoproterenol(n);
if strcmp(checkFigure, '2b')
p(17) = IBMX(n);
plot_cAMP = true;
end
case {'2c','2d'}
y0(16) = cAMP_holding(n); %cAMP
y0(15) = cAMP_holding(n)*4;
p(32) = PKI(n);
case {'2f', '2gh'}
p(1) = isoproterenol(n);
case {'3a', '3b'}
plot_cAMP = true;
p(17) = IBMX(n);
p(1) = isoproterenol(n);
if n > numloops0 && n <= numloops0*2
p(14) = p_orig(14)*expr_mult(n); % AC_tot
elseif n > numloops0*2 && n <= numloops0*3
y0(4) = y0_orig(4)*expr_mult(n); % b1AR_tot
elseif n > numloops0*3 && n <= numloops0*4
p(27) = p_orig(27)/expr_mult(n); % Kgsa
end
case {'4a','4b'}
plot_cAMP = true;
p(1) = isoproterenol(n);
p(17) = IBMX(n);
p(12) = p_orig(12)*(0.05*CTX(n) + ~CTX(n)); % p_12 = k_hyd
case '5'
plot_cAMP = true;
p(1) = isoproterenol(n);
if n > numloops0 && n <= numloops0*2
p(30) = p_orig(30)*0.05; % PKAItot
elseif n > 2*numloops0 && n <= numloops0*3
p(32) = p_orig(32)*0.05; % PKItot
end
case '6a'
plot_cAMP = true;
p(1) = isoproterenol(n);
if n == 2
fig6p = 'LCC';
y0(19) = 0;
elseif n == 3
fig6p = 'PLB';
y0(23) = 0;
y0(24) = 0;
end
case '6b'
p(1) = isoproterenol(n);
PLBKO = vPLBKO(n);
if vPLBKO(n) && false
p(38) = p_orig(38)*0.00001; % PLBtot
y0(19) = y0_orig(19)*0.00001; % PLBp
end
end
options = odeset('Mass',M,'RelTol',1e-6,'MaxStep',5e-3,'Stats','off');
% -------------------- Run simulation --------------------
stp = struct();
stp.protocol = protocol;
stp.p = p;
stp.fig6 = fig6p;
stp.PLBKO = PLBKO;
stp.fracPLBo = fracPLBo;
stp.checkFigure = checkFigure;
[t_temp,y] = ode15s(@f,tspan,y0,options,stp);
[~,datOut] = f(t_temp, y, stp);
of_2fgh(n) = false;
if strcmp(checkFigure, '2f') || strcmp(checkFigure, '2gh')
of_2fgh(n) = true;
end
if (n == 1 && of_2fgh(n)) || ~of_2fgh(n)
data_t1(n,:) = t_temp;
data_y1(n,:,:) = y;
data_I1(n,:,:) = [datOut.I_Ca, datOut.I_app];
end
switch checkFigure
case {'2f', '2gh'}
if n == 2
data_t2(n,:) = t_temp;
data_y2(n,:,:) = y;
data_I2(n,:,:) = [datOut.I_Ca, datOut.I_app];
end
case {'2c','2d'}
if n <= numloops0
data_cAMP_ss0(n,:) = [y(end,15), y(end,16)];
data_PKA0(n,:) = y(end,17);
else
data_cAMP_ss1(n,:) = [y(end,15), y(end,16)];
data_PKA1(n,:) = y(end,17);
end
case '2e'
data_fracPLBp = [data_fracPLBp, datOut.fracPLBp(end)];
case '2b'
if n <= numloops0
data_cAMP_ss_IBMX0(n) = y(end,15); %[1:int:end]
else
data_cAMP_ss_IBMX1(n) = y(end,15);
end
case {'3a', '3b'}
if n <= numloops0
data_cAMP_normal = [data_cAMP_normal, y(end,15)];
elseif n <= 2*numloops0 && n > numloops0
data_cAMP_ACmult = [data_cAMP_ACmult, y(end,15)];
elseif n <= 3*numloops0 && n > 2*numloops0
data_cAMP_B1ARmult = [data_cAMP_B1ARmult, y(end,15)];
else
data_cAMP_KGsadiv = [data_cAMP_KGsadiv, y(end,15)];
end
case {'4a','4b'}
data_fig4 = [data_fig4, y(end,15)]; % cAMP_tot
data_Gsa_gtp_AC(n) = datOut.Gsa_gtp_AC(end);
case '5'
if n <= numloops0
data_cAMP_Control = [data_cAMP_Control, y(end,15)];
elseif n <= 2*numloops0 && n > numloops0
data_cAMP_PKA_KO = [data_cAMP_PKA_KO, y(end,15)];
elseif n <= 3*numloops0 && n > 2*numloops0
data_cAMP_PKI_KO = [data_cAMP_PKI_KO, y(end,15)];
end
case '6a'
data_CaTransient = [data_CaTransient; y(:,42)'];
case '6b'
iSS = find(t_temp>tspan(end)-1,1);
data_CaTransient_PLBKO = [data_CaTransient_PLBKO, max(y(iSS:end,42))];
end
% plot cAMP timecourse in TEMP folder
if true
loop_plot = false;
ip = zeros(numloops,1);
legnames = cell(numloops,1);
pm = cAMP_mult;
ycaption = 'conc [pmol/mg]';
if plot_cAMP
ip(n) = 15;
legnames{n} = 'cAMP_tot';
loop_plot = true;
elseif (strcmp(checkFigure, '2c') || strcmp(checkFigure, '2d'))
ip(n,1) = 15;
legnames{n} = 'cAMP_tot';
loop_plot = true;
elseif strcmp(checkFigure, '2e')
ip(n,1) = 19;
legnames{n} = {'PLBp'};
pm = 1/p_orig(38);
ycaption = 'conc uM';
loop_plot = true;
elseif strcmp(checkFigure, '6b')
ip(n,1) = 42;
legnames{n} = {'calcium_i mM'};
ycaption = 'conc uM';
pm = 1;
loop_plot = true;
end
if loop_plot
figure,
for j=1:length(ip(n))
plot(t_temp, y(:,ip(n,j))*pm); hold on;
end
hold off;
legend(legnames{n});
xlabel('time')
ylabel(ycaption)
title(sprintf('loop %d',n));
saveas(gcf,[temp_path num2str(n) '.jpg']);
end
end
end
I_Ca1 = data_I1(1,:,1);
I_app1 = data_I1(1,:,2);
Cai1 = data_y1(1,:,42);
Vm1 = data_y1(1,:,43);
if of_2fgh
I_Ca2 = data_I2(2,:,1);
I_app2 = data_I2(2,:,2);
Cai2 = data_y2(2,:,42);
Vm2 = data_y2(2,:,43);
end
if exist('data_t1','var')
% translate time for experiments under protocol 1 or 2
data_t1 = data_t1 - 59;
end
if exist('data_t2','var')
% translate time for experiments under protocol 1 or 2
data_t2 = data_t2 - 59;
end
co = [0 0 1
1 0 0
0 1 0];
set(groot,'defaultAxesColorOrder',co, 'DefaultLineLineWidth', 4, 'DefaultLineMarkerSize',12)
figure('DefaultAxesFontSize',18)
switch checkFigure
case '2a'
dat = csvread([extracted_figure_path checkFigure '.csv'],1,0); % skip header
plot(dat(:,1),dat(:,2),'bx'); hold on;
plot(data_t1/60, data_y1(1,:,15)*cAMP_mult,'b');
xlabel ('Time (min)')
ylabel('cAMP_{tot} (pmol/mg)');
% legend('primary data','reproduction');
case '2b'
datwithIBMX = csvread([extracted_figure_path '2b_withIBMX.csv'],1,0); % skip header
datwithoutIBMX = csvread([extracted_figure_path '2b_withoutIBMX.csv'],1,0); % skip header
semilogx(isoproterenol(1:numloops0), data_cAMP_ss_IBMX0(1:numloops0)*cAMP_mult, 'b-'); hold on;
semilogx(isoproterenol(numloops0+1:end), data_cAMP_ss_IBMX1(numloops0+1:end)*cAMP_mult, 'r-')
semilogx(datwithoutIBMX(:,1),datwithoutIBMX(:,2),'bx'); hold on;
semilogx(datwithIBMX(:,1),datwithIBMX(:,2),'rx'); hold on;
legend('no IBMX','with IBMX','Location','east');
xlim([10e-3 10e1]);
xlabel('Isoproterenol (\muM)')
ylabel('cAMP_{tot} (pmol/mg)')
case {'2c','2d'}
denom = 2*p_orig(30); % PKAI_tot
F1 = data_PKA0(1:numloops0)/denom;
F2 = data_PKA1(numloops0+1:end)/denom;
datSolid = csvread([extracted_figure_path checkFigure '_solid.csv'],1,0); % skip header
datDashed = csvread([extracted_figure_path checkFigure '_dashed.csv'],1,0);
if strcmp(checkFigure,'2c')
semilogx(data_cAMP_ss0(1:numloops0,1),F1,'b'); hold on;
semilogx(data_cAMP_ss1(numloops0+1:end,1),F2,'r');
semilogx(datSolid(:,1),datSolid(:,2),'bx'); hold on;
semilogx(datDashed(:,1),datDashed(:,2),'rx'); hold on;
hold off;
legend('PKACI', 'PKACI+PKI','Location','northwest');
xlabel('cAMP_{tot} (\muM)')
ylabel('PKA activation')
else
semilogx(data_cAMP_ss0(1:numloops0,1),log10(F1./(1-F1)),'b'); hold on;
semilogx(data_cAMP_ss1(numloops0+1:end,1),log10(F2./(1-F2)),'r'); hold on;
semilogx(datSolid(:,1),datSolid(:,2),'bx'); hold on;
semilogx(datDashed(:,1),datDashed(:,2),'rx'); hold on;
hold off;
legend('PKACI', 'PKACI+PKI','Location','northwest');
xlabel('cAMP_{tot} (\muM)')
ylabel('log(F/(1-F))')
xlim([0.01 1]);
end
case '2e'
dat = csvread([extracted_figure_path checkFigure '.csv'],1,0); % skip header
semilogx(dat(:,1),dat(:,2),'bx'); hold on;
semilogx(isoproterenol, data_fracPLBp,'b'); %hold on;
xlabel('Isoproterenol (\muM)');
ylabel('frac PLBp');
% legend('primary data','reproduction','Location','east');
case '2f'
datSolid = csvread([extracted_figure_path checkFigure '_solid.csv'],1,0); % skip header
datDotted = csvread([extracted_figure_path checkFigure '_dotted.csv'],1,0); % skip header
queryTime = tspan(end)-1;
iSS = [find(tspan == queryTime), find(tspan == queryTime+0.3), length(tspan)];
plot(data_t1(1,iSS(1):iSS(2)),I_Ca1(iSS(1):iSS(2))*Cm*1e3,'b'); hold on;
plot(data_t2(2,iSS(1):iSS(2)),I_Ca2(iSS(1):iSS(2))*Cm*1e3,'r');
plot(datSolid(:,1),datSolid(:,2),'bx'); hold on;
plot(datDotted(:,1),datDotted(:,2),'rx'); hold on;
legend('no Iso','1 \muM Iso', 'Location','southeast')
ylabel('I_{Ca} (nA)')
xlabel('Time (s)')
case '2gh' % query Figs 2g-h
datSolid = csvread([extracted_figure_path '2g_solid.csv'],1,0); % skip header
datDashed = csvread([extracted_figure_path '2g_dashed.csv'],1,0);
queryTime = tspan(end)-1;
iSS = [find(tspan == queryTime), find(tspan == queryTime+0.3), length(tspan)];
plot(data_t1(1,iSS(1):iSS(3)),Cai1(iSS(1):iSS(3))*1e3,'b'); hold on;
plot(data_t2(2,iSS(1):iSS(3)),Cai2(iSS(1):iSS(3))*1e3,'r'); hold on;
plot(datSolid(:,1),datSolid(:,2),'bx'); hold on;
plot(datDashed(:,1),datDashed(:,2),'rx'); hold on;
ylabel('[Ca_i] (\muM)')
xlabel('Time (s)')
legend('no Iso', '1 \muM Iso','Location','northeast');
title('2g')
saveas(gcf,'figures/F2g.jpg');
% subplot(1,2,2)
set(groot,'defaultAxesColorOrder',co, 'DefaultLineLineWidth', 4, 'DefaultLineMarkerSize',12)
figure('DefaultAxesFontSize',18)
datSolid = csvread([extracted_figure_path '2h_solid.csv'],1,0); % skip header
datDashed = csvread([extracted_figure_path '2h_dashed.csv'],1,0);
plot(data_t1(iSS(1):iSS(2)), Vm1(iSS(1):iSS(2)),'b'); hold on;
plot(data_t2(2,iSS(1):iSS(2)), Vm2(iSS(1):iSS(2)),'r'); hold on;
plot(datSolid(:,1),datSolid(:,2),'bx'); hold on;
plot(datDashed(:,1),datDashed(:,2),'rx'); hold on;
legend('no Iso', '1 \muM Iso','Location','northeast');
xlabel('Time (s)')
ylabel('Membrane Voltage (mV)')
ylim([-100 20])
title('2h')
saveas(gcf,'figures/F2h.jpg');
case {'3a', '3b'}
datSolid = csvread([extracted_figure_path checkFigure '_solid.csv'],1,0); % skip header
datDashed = csvread([extracted_figure_path checkFigure '_dashed.csv'],1,0); % skip header
if strcmp(checkFigure, '3a')
datmidDashed = csvread([extracted_figure_path checkFigure '_midDashed.csv'],1,0); % skip header
datDotted = csvread([extracted_figure_path checkFigure '_dotted.csv'],1,0); % skip header
end
semilogx(isoproterenol(1:numloops0), data_cAMP_normal*cAMP_mult, 'b-'); hold on;
semilogx(isoproterenol(numloops0+1:numloops0*2), data_cAMP_ACmult*cAMP_mult, 'r-'); hold on;
semilogx(isoproterenol(numloops0*2+1:numloops0*3), data_cAMP_B1ARmult*cAMP_mult, 'g-'); hold on;
semilogx(isoproterenol(numloops0*3+1:end), data_cAMP_KGsadiv*cAMP_mult, 'm-'); hold on;
semilogx(datSolid(:,1),datSolid(:,2),'bx'); hold on;
semilogx(datDashed(:,1),datDashed(:,2),'rx'); hold on;
if strcmp(checkFigure, '3a')
semilogx(datDotted(:,1),datDotted(:,2),'gx'); hold on;
end
if strcmp(checkFigure, '3a')
semilogx(datmidDashed(:,1),datmidDashed(:,2),'mx'); hold on;
end
if strcmp(checkFigure, '3a')
xlim([10e-4,1]);
ylim([0 650]);
legend(...
{'Control', ...
sprintf('[AC_{tot}] * %g', expr_mult0(2)), ...
sprintf('[\\beta1AR_{tot}] * %g', expr_mult0(3)), ...
sprintf('K_{Gs\\alpha} / %g', expr_mult0(4)),...
},...
'Location','northwest','FontSize',14);
else
legend(...
{'Control', ...
sprintf('[AC_{tot}] * %g', expr_mult0(2)), ...
sprintf('[\\beta1AR_{tot}] * %g', expr_mult0(3)), ...
sprintf('KGsa / %g', expr_mult0(4))},...
'Location','northwest','FontSize',14);
end
xlabel('Isoproterenol (\muM)')
ylabel('cAMP_{tot} (pmol/mg)')
case '4a'
dat = csvread([extracted_figure_path checkFigure '.csv'],1,0); % skip header
denom = data_fig4(1); %normalise by data at Iso=1uM, CTX=0
data_fig4 = data_fig4/denom;
bar([dat';data_fig4]',1)
ylabel('AC activity') % cAMP_{tot} (\muM)
set(gca, 'XTickLabel',{'Iso','CTX','Iso+CTX'});
legend('primary data','reproduction','Location','northwest');
case '4b'
datSolid = csvread([extracted_figure_path checkFigure '_solid.csv'],1,0); % skip header
datDotted = csvread([extracted_figure_path checkFigure '_dotted.csv'],1,0); % skip header
denom = data_fig4(end); % normalise by data at 10uM ISO+CTX
y = data_fig4/denom;
x = data_Gsa_gtp_AC/p_orig(14); % p not available after for loop
plot(x(1:2),y(1:2),'b-'); hold on;
plot(x(3:4),y(3:4),'r-');
plot(datSolid(:,1),datSolid(:,2),'bx'); hold on;
plot(datDotted(:,1),datDotted(:,2),'rx'); hold on;
legend('Control','CTX','Location','northwest');
xlabel('Fractional AC activation'); % GsaGTP_{AC}/AC_{tot}
ylabel('f_{cAMP}') % cAMP_{tot}/cAMP_{tot,max}
for j=1:numloops
protostr = sprintf('%g \\muM', isoproterenol(j));
text(x(j)+0.01,y(j)-0.03,protostr,'fontsize',16)
end
case '5'
dat = csvread([extracted_figure_path checkFigure '_Control.csv'],1,0); % skip header
datPKI = csvread([extracted_figure_path checkFigure '_PKI.csv'],1,0);
datPKA = csvread([extracted_figure_path checkFigure '_PKA.csv'],1,0);
semilogx(isoproterenol(1:numloops0), data_cAMP_Control*cAMP_mult, 'b-'); hold on;
semilogx(isoproterenol(1:numloops0), data_cAMP_PKA_KO*cAMP_mult, 'r-'); hold on;
semilogx(isoproterenol(1:numloops0), data_cAMP_PKI_KO*cAMP_mult, 'g-');
semilogx(dat(:,1),dat(:,2),'bx'); hold on;
semilogx(datPKA(:,1),datPKA(:,2),'rx'); hold on;
semilogx(datPKI(:,1),datPKI(:,2),'gx'); hold on;
legend('Control',...
'PKA-KO', ...
'PKI-KO', 'Location','east');
xlabel('Isoproterenol (\muM)')
ylabel('cAMP_{tot} (pmol/mg)')
case '6a'
dat = csvread([extracted_figure_path checkFigure '_both.csv'],1,0); % skip header
datLCC = csvread([extracted_figure_path checkFigure '_LCC.csv'],1,0);
datPLB = csvread([extracted_figure_path checkFigure '_PLB.csv'],1,0);
queryTime = tspan(end)-1;
iSS = [find(tspan == queryTime), length(tspan)];
plot(data_t1(1,iSS(1):iSS(2)), data_CaTransient(1,(iSS(1):iSS(2)))*1e3,'b-'); hold on; % convert mM to uM
plot(data_t1(1,iSS(1):iSS(2)), data_CaTransient(2,(iSS(1):iSS(2)))*1e3,'r-'); hold on;
plot(data_t1(1,iSS(1):iSS(2)), data_CaTransient(3,(iSS(1):iSS(2)))*1e3,'g-'); hold on;
plot(dat(:,1),dat(:,2),'bx'); hold on;
plot(datLCC(:,1),datLCC(:,2),'rx'); hold on;
plot(datPLB(:,1),datPLB(:,2),'gx'); hold on;
legend('PLB + LCC', 'LCC','PLB','Location','northeast');
xlabel('Time (s)')
ylabel('Calcium (\muM)');
hold off;
case '6b'
dat = csvread([extracted_figure_path checkFigure '.csv'],1,0); % skip header
d = [data_CaTransient_PLBKO*1e3];
bar([dat';d]', 1)
ylabel('Calcium (\muM)');
set(gca, 'XTickLabel',{'None','Iso','PLB KO','Iso+PLB KO'});
legend('primary data','reproduction','Location','northwest');
otherwise
disp('figure not specified properly');
end
if ~strcmp(checkFigure, '2gh')
title(checkFigure);
saveas(gcf,['figures/F' checkFigure '.jpg']);
end
end
