Cardiac models of circadian rhythms in early afterdepolarizations & arrhythmias (Diekman & Wei 2021)

 Download zip file 
Help downloading and running models
Accession:266974
We fit a simplified Luo-Rudy model to voltage-clamp data showing a circadian rhythm in L-type calcium conductance. Simulations of the model (single-cell and 2-D spatial versions) suggest that circadian rhythms in early afterdepolarizations may contribute to daily rhythms in cardiac arrhythmias and sudden cardiac death.
Reference:
1 . Diekman CO, Wei N (2021) Circadian Rhythms of Early Afterdepolarizations and Ventricular Arrhythmias in a Cardiomyocyte Model. Biophys J 120:319-333 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Cardiac ventricular cell;
Channel(s): I Calcium; I K;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: MATLAB;
Model Concept(s): Circadian Rhythms; Action Potentials; Cardiac pacemaking; Bifurcation;
Implementer(s): Diekman, Casey O. [casey.o.diekman at njit.edu];
Search NeuronDB for information about:  I K; I Calcium;
clear all; close('all')

set(0,'DefaultAxesFontSize',24)

StimProtocol=2;

% This script is based on Peter Hammer's code "Spiral waves in monodomain
% reaction-diffusion-model" provided on the MATLAB Central File Exchange:
% https://uk.mathworks.com/matlabcentral/fileexchange/22492-spiral-waves-in-monodomain-reaction-diffusion-model

% Matlab implementation a monodomain reaction-diffusion model in 2-D. 
% The model equations are a variant of the Fitzhugh-Nagumo equations
% modified to simulate the cardiac action potential. The progression
% of the two normalized state variables, membrane voltage (v) and recovery 
% (r), is computed across a 128 x 128 spatial domain and across time. This 
% function simulates spiral waves, which are hypothesized to underlie 
% reentrant tachycardia. The spiral waves can be initiated by two different 
% cardiac pacing methods: 
%
% (1) two-point stimulation where a point stimulus is delivered in the 
% center of the domain followed by another point stimulus on the partially
% refractory wake of the first wave of excitation.
%
% (2) cross-field stimulation where a stimulus is applied to the left 
% domain boundary causing a plane wave. As this wave travels across the
% domain, a second stimulus is applied to the bottom boundary of the domain.
%
% This function accepts only one input argument, StimProtocol, which can 
% be either the numerical values '1' (for two-point stimulation) or '2' 
% (for cross-field stimulation). As the simulation runs, the activation 
% state of the individual units comprising the domain is mapped to color 
% and plotted in a figure window. A count of time steps is displayed at 
% the top of the plot along with the values of v and r for the upper left 
% element of the domain. 
% 
% Model equations are solved using a finite difference method for spatial 
% derivatives and explicit Euler integration for time derivatives. Newman
% boundary conditions are assumed. Model parameters are taken from two 
% journal articles: [1] Rogers JM et al. "A collocation-Galerkin finite 
% element model of cardiac action potential propagation". IEEE TBME;41:
% 743-757,1994. [2] Pertsov AM et al. "Spiral waves of excitation underlie 
% reentrant activity in isolated cardiac muscle". Circulation Research;72:
% 631-650, 1993. 
%
% Following the simulated spiral waves, a movie (AVI) is generated, and 
% the user is given the option to save the movie to disk. One simulation 
% takes about 160 seconds on a 2.33 GHz Intel Dual Core 64-bit laptop 
% (3.3 GB RAM).
%
% This function was written by Peter E. Hammer (hammer@usa.com) in 
% October 2006 as part of an academic project associated with PhD studies 
% in biomedical engineering. Have fun!

ncols=128;                               % Number of columns in domain
nrows=128;                               % Number of rows in domain
dur=60000;                               % Number of time steps
h=0.01;                                 % Grid size
h2=h^2;
dt=0.1;                                     % Time step                                
Iex=500;                                % Amplitude of external current
mu=1.0;                                  % Anisotropy
Gx=0.000025; Gy=Gx/mu;                  % Conductances               
v=-80*ones(nrows,ncols);                % Initialize voltage array
fgate=ones(nrows,ncols);                % Initialize gatingg variable array
xgate=zeros(nrows,ncols);               % Initialize gatingg variable array

% we save the potential at several points for plotting
% this is the index to plot the potential at (center of domain)
rowInd=round(nrows/2);
colInd=round(ncols/2);

% Parameters for Diekman and Wei model
Iapp=0;
Eca=60;
Ek=-80;
%g_ca=0.15;
g_ca=0.3;
g_k=0.1;
%g_k=0.05;

gk_mat=g_k*ones(nrows,ncols);
gk_mat(:,1:102)=0.05;                 % heterogeniety in potassium channel conductance

dhalf=7.3; dslope=8.6; fhalf=13.3; fslope=11.9;
xhalf = 40; xslope = 5;

tauf=80;
taux=300;

% Set initial stim current and pattern
iex=zeros(nrows,ncols);
if StimProtocol==1
    iex(62:67,62:67)=Iex;
else
    iex(:,1)=Iex;
end

clim=[-80 60];
% Setup image
ih=imagesc(v,clim); %set(ih,'cdatamapping','direct')
colormap(hot); axis image off; th=title('');
colorbar

set(gcf, 'PaperPosition', [0.25 2.5 8.0 6.0]);


n=0;                 % Counter for time loop
k=0;                 % Counter for movie frames
done=0;              % Flag for while loop

n1e=20;              % Step at which to end 1st stimulus
switch StimProtocol
    case 1           % Two-point stimulation
        n2b=3800;    % Step at which to begin 2nd stimulus
        n2e=3900;    % Step at which to end 2nd stimulus
    case 2           % Cross-field stimulation
        n2b=8100;     % Step at which to begin 2nd stimulus
        n2e=8130;     % Step at which to end 2nd stimulus
end

n3b=10000;
n3e=10030;

while ~done          % Time loop
    
    if n == n1e      % End 1st stimulus
        iex=zeros(nrows,ncols);
    end
    
    if n == n2b      % Begin 2nd stimulus
        switch StimProtocol
            case 1
                iex(62:67,49:54)=Iex;
            case 2
                iex(end,:)=Iex;
        end
    end
    
    if n == n2e      % End 2nd stimulus
        iex=zeros(nrows,ncols);
    end
    
    % Create padded v matrix to incorporate Newman boundary conditions 
    vv=[[0 v(2,:) 0];[v(:,2) v v(:,end-1)];[0 v(end-1,:) 0]];
    
    % Update v
    vxx=(vv(2:end-1,1:end-2) + vv(2:end-1,3:end) -2*v)/h2; 
    vyy=(vv(1:end-2,2:end-1) + vv(3:end,2:end-1) -2*v)/h2;
        
    dinf = 1./(1+exp(-(v+dhalf)/dslope));
    finf = 1./(1+exp((v+fhalf)/fslope));
    
    Ica = g_ca*dinf.*fgate.*(v-Eca);
        
    xinf = 1./(1+exp(-(v+xhalf)/xslope));
    %Ik = g_k*xgate.*(v-Ek);
    Ik = gk_mat.*xgate.*(v-Ek);
    
    dvdt = Iapp - Ica - Ik + iex + Gx*vxx + Gy*vyy;

    v_new = v + dvdt*dt;
    
    dfdt=(finf-fgate)/tauf;
    fgate=fgate+dfdt*dt;
    dxdt=(xinf-xgate)/taux;
    xgate=xgate+dxdt*dt;
    v=v_new; clear v_new
    
    VsaveMiddleCol(n+1,:)=v(:,64);
    VsaveMiddleRow(n+1,:)=v(64,:);
    
    % Update image and text 
    set(ih,'cdata',v);
    set(th,'string',sprintf('%d  %0.2f   %0.2f    %0.2f',n,v(rowInd,colInd),fgate(rowInd,colInd),xgate(rowInd,colInd)))
    drawnow
        
    n=n+1;
    done=(n > dur);

end


%% make plot similar to Figure 7D

t=(0:(n-1))*dt;

figure(2)
hold on
plot(t,VsaveMiddleRow(:,1),'b','linewidth',3)
plot(t,VsaveMiddleRow(:,64),'r','linewidth',3)
plot(t,VsaveMiddleRow(:,128),'k','linewidth',3)
ylim([-90 60])
set(gca,'XTick',0:1000:5000,'YTick',[-80:40:80])
xlabel('$t$ (ms)','interpreter','latex')
ylabel('$V$ (mV)','interpreter','latex')




Loading data, please wait...