-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathFigure6.m
112 lines (96 loc) · 4.76 KB
/
Figure6.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
%% Figure 6:AP clamp simulations illustrate changes in K+ currents with
%% alterations in AP shape.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--- "Slow delayed rectifier current protects ventricular myocytes from
% arrhythmic dynamics across multiple species: a computational study" ---%
% By: Varshneya,Devenyi,Sobie
% For questions, please contact Dr.Eric A Sobie -> [email protected]
% or put in a pull request or open an issue on the github repository:
% https://github.com/meeravarshneya1234/IKs_stabilizes_APs.git.
%--- Note:
% Results displayed in manuscript were run using MATLAB 2016a on a 64bit
% Intel Processor. For exact replication of figures it is best to use these
% settings.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------------------------------------------------------------
%% Figure 6
%--- Description of Figure:
% AP Clamp simulations performed in the O'Hara model
%---: Functions required to run this script :---%
% main_program.m - runs single AP simulation
% main_APClamp.m - runs the AP clamp simulation using parfor loop
% plotting_APClamp.m - plots the waveforms from the main_APClamp.m function
%--------------------------------------------------------------------------
%%
%---- Run under baseline conditions to get input voltage and time for AP
% Clamp ----%
settings.model_name = 'Ohara';
settings.celltype = 'endo';
settings.PCL = 1000 ; % Interval bewteen stimuli,[ms]
settings.stim_delay = 100 ; % Time the first stimulus, [ms]
settings.stim_dur = 2 ; % Stimulus duration
settings.stim_amp = 32.2; % Stimulus amplitude
settings.nBeats = 1 ; % Number of beats to simulate
settings.numbertokeep = 1;% Determine how many beats to keep. 1 = last beat, 2 = last two beats
settings.steady_state = 1;
% Ks and Kr scaling must have the same length vector.
settings.Ks_scale = 1; % perturb IKs, set to 1 for baseline
settings.Kr_scale = 1; % perturb IKr, set to 1 for baseline
settings.Ca_scale =1; % perturb ICaL, set to 1 for baseline
base_datatable= main_program(settings);
statevars_input = base_datatable.states{1,1}; % used as input for the AP clamp
t_input = base_datatable.times{1,1}; % used as input for the AP clamp
%---- Run AP Clamp on a linear scale ----%
apclamp_settings.model_name = 'Ohara';
apclamp_settings.celltype = 'endo';
[~,Vind] = ICs(apclamp_settings.model_name,1,1000);
statevar_i_original = statevars_input(1:end ~= Vind); % all state variables besides V
apclamp_settings.statevar_i_original = cellfun(@(v) v(end), statevar_i_original);
apclamp_settings.volts = statevars_input{1,Vind};
apclamp_settings.times = t_input;
apclamp_settings.repol_change = [0.5,1,1.5,2,2.5]; % scaling
apclamp_settings.numberofAPs = 100;
apclamp_settings.numberkeep = 1;
apclamp_settings.PCL = 1000;
apclamp_settings.vals_graph = [0.5,1,1.5,2,2.5];
APClamp_data = main_APClamp(apclamp_settings); % run simulation
indxs = find(ismember(round(apclamp_settings.repol_change,2),round(apclamp_settings.vals_graph,2)));
data_to_graph = APClamp_data(:,indxs);
%% Plot Figure 6A & 6B
plotting_APClamp(apclamp_settings,data_to_graph) % plot results
%---- Run AP Clamp on a log scale ----%
apclamp_settings.model_name = 'Ohara';
apclamp_settings.celltype = 'endo';
[~,Vind] = ICs(apclamp_settings.model_name,1,1000);
statevar_i_original = statevars_input(1:end ~= Vind); % all state variables besides V
apclamp_settings.statevar_i_original = cellfun(@(v) v(end), statevar_i_original);
apclamp_settings.volts = statevars_input{1,Vind};
apclamp_settings.times = t_input;
logscalefactors = linspace(log(1/3),log(3),11) ;
apclamp_settings.repol_change = exp(logscalefactors) ; % scaling
apclamp_settings.numberofAPs = 100;
apclamp_settings.numberkeep = 1;
apclamp_settings.PCL = 1000;
log_APClamp_data = main_APClamp(apclamp_settings);
%% Plot Figure 6C
BL = find(round(apclamp_settings.repol_change,2) ==1.00);
AKr = cell2mat(arrayfun(@(x) x.Area_Kr(1), log_APClamp_data,'UniformOutput', 0));
AKs = cell2mat(arrayfun(@(x) x.Area_Ks(1), log_APClamp_data,'UniformOutput', 0));
Volts = cell2mat(arrayfun(@(x) x.Voltage(:,1), log_APClamp_data,'UniformOutput', 0));
times = cell2mat(arrayfun(@(x) x.time(:,1), log_APClamp_data,'UniformOutput', 0));
IKr = cell2mat(arrayfun(@(x) x.IKr(:,1), log_APClamp_data,'UniformOutput', 0));
IKs = cell2mat(arrayfun(@(x) x.IKs(:,1), log_APClamp_data,'UniformOutput', 0));
figure
loglog(apclamp_settings.repol_change,AKr/AKr(BL),'k-x', 'linewidth',2)
hold on
loglog(apclamp_settings.repol_change,AKs/AKs(BL),'k-o','linewidth',2)
xlabel('APD Factor')
ylabel('norm AUC')
legend('IKr','IKs')
set(gca,'FontSize',12,'FontWeight','bold')
title('Figure 6C')
figure
bar([AKr(end)/AKr(BL),AKs(end)/AKs(BL)],0.5)
ylabel('norm current (A/F)')
set(gca,'XTickLabels',{'IKr','IKs'})
title('Figure 6C - inset')