Circadian rhythms of early afterdepolarizations and ventricular arrhythmias in a cardiomyocyte model

Sudden cardiac arrest is a malfunction of the heart’s electrical system, typically caused by ventricular arrhythmias, that can lead to sudden cardiac death (SCD) within minutes. Epidemiological studies have shown that SCD and ventricular arrhythmias are more likely to occur in the morning than in the evening, and laboratory studies indicate that these daily rhythms in adverse cardiovascular events are at least partially under the control of the endogenous circadian timekeeping system. However, the biophysical mechanisms linking molecular circadian clocks to cardiac arrhythmogenesis are not fully understood. Recent experiments have shown that L-type calcium channels exhibit circadian rhythms in both expression and function in guinea pig ventricular cardiomyocytes. We developed an electrophysiological model of these cells to simulate the effect of circadian variation in L-type calcium conductance. We found that there is a circadian pattern in the occurrence of early afterdepolarizations (EADs), which are abnormal depolarizations during the repolarization phase of a cardiac action potential that can trigger fatal ventricular arrhythmias. Specifically, the model produces EADs in the morning but not at other times of day. We show that the model exhibits a codimension-2 Takens-Bogdanov bifurcation that serves as an organizing center for different types of EAD dynamics. We also simulated a 2-D spatial version of this model across a circadian cycle. We found that there is a circadian pattern in the breakup of spiral waves, which represents ventricular fibrillation in cardiac tissue. Specifically, the model produces spiral wave breakup in the morning but not in the evening. Our study is the first to propose a link between circadian rhythms and EAD formation and suggests that the efficacy of drugs targeting EAD-mediated arrhythmias may depend on the time of day that they are administered. Significance Statement Why are life-threatening cardiac arrhythmias more likely to occur in the morning than in the evening? The electrical properties of the heart exhibit daily rhythms due to molecular circadian clocks within cardiomyocytes. Our computational model of ventricular myocytes shows that clock-controlled expression of a voltage-gated calcium ion channel leads to early afterdepolarizations (EADs) at certain times of the day. EADs, in which the membrane potential of a cardiomyocyte depolarizes a second time before fully repolarizing, can trigger arrhythmias. To our knowledge, this is the first study linking the circadian clock to EAD formation. Our results suggest that the efficacy of anti-arrhythmic medications targeting this calcium ion channel may depend on the time of day the drug is taken.


Introduction
Sudden cardiac arrest (SCA) is the most common single cause of natural death in the United States [62]. Distinct from a heart attack, SCA occurs when the electrical system of the heart malfunctions, often without prior symptoms. It is usually caused by arrhythmias such as ventricular tachycardia and ventricular fibrillation. These abnormally fast and irregular heartbeats do not pump blood properly and can cause sudden cardiac death (SCD) within minutes if emergency treatment is not begun immediately [61].
The risk of sudden cardiac arrest is not constant throughout the day. SCD is more likely to occur in the morning than in the evening [47,72]. Ventricular tachyarrhythmias also exhibit a diurnal rhythm with a peak in the morning [32,56]. The biophysical mechanisms underlying these daily rhythms in adverse cardiovascular events are not fully understood. The master circadian (∼24-hour) pacemaker in the hypothalamus, the suprachiasmatic nucleus (SCN), influences a variety a cardiovascular phenomenon by coordinating daily rhythms in the release of hormones and other circulating molecules. Recently, it has been demonstrated that circadian clocks within heart muscle cells (cardiomyocytes) also regulate rhythms in cardiac electrophysiology [3].
These intracellular circadian clocks are comprised of transcriptional and translational feedback loops that lead to ∼24-hour rhythms in gene expression. In mice, cardiac ion channel expression and myocardial repolarization are under the control of a clock-dependent oscillator that regulates potassium channel-interacting protein 2 (KChIP2), a subunit required for generating the transient outward potassium current I to [29]. Reduced I to amplitude has arrhythmogenic consequences, perhaps due to lengthened QT (repolarization and depolarization) intervals, and may contribute to sudden death in the early stages of human heart failure [25]. The effect of circadian variation in potassium current on action potential (AP) duration and QT interval has been studied using mathematical models of murine, guinea pig, and human myocytes [21,60]. Recent experiments in guinea pig myocytes have shown that L-type calcium current (I CaL ) is under circadian control as well, possibly through the PI3K-Akt signaling pathway [9]. How QT interval is affected by circadian oscillations in the concentration of sodium, potassium, and calcium ions in plasma was also studied in biophysically detailed models of human left ventricular cardiomyocytes [20].
In addition to lengthened QT intervals, the presence of early afterdepolarizations (EADs) is also associated with the development of ventricular arrhythmias [71]. EADs are voltage deflections that occur before full repolarization of the membrane potential during an AP. Extensive modeling of EADs has been performed to understand the ionic and dynamical mechanisms involved in the generation of EADs in isolated cells and their spatial propagation in cardiac tissue [33,73,77,68,34,70]. At a basic level, EADs result from reduced repolarization reserve due to reduced outward potassium currents or elevated inward calcium currents [53]. Thus, circadian variation in these currents could render myocytes more vulnerable to EADs at certain times of day, and play a role in the observed circadian profile of ventricular arrhythmias and SCDs.
In this paper, we use biophysical modeling and dynamical systems analysis to study how circadian variation in ionic conductances affects EAD generation. First, we fit a conductance-based model to published electrophysiological data from guinea pig ventricular myocytes at two circadian time points. We then perform simulations of single-cell and 2-D spatial domain versions of the model across a circadian cycle. In the single-cell model, we find that EADs occur in the morning but not at other times of day. In the spatial model, we observe that spiral wave breakup, a phenomenon associated with ventricular arrhythmias in cardiac tissue, occurs in the morning but not in the evening. We also show that the single-cell model exhibits a codimension-2 Takens-Bogdanov bifurcation, which can serve as an organizing center for the different types of EAD dynamics that have been observed. To the best of our knowledge, this work is the first to consider connections between the circadian clock and EADs.

Single-cell model
We used recently published voltage-clamp recordings from guinea pig ventricular myocytes to modify the Sato et al. [57] minimal model of cardiac action potential generation. This conductance-based model describes the dynamics of the membrane potential V using the Hodgkin-Huxley modeling formalism and is a three-dimensional system of ordinary differential equations (ODEs): The model includes an inward L-type calcium current (I CaL ), an outward potassium current (I K ), and an externally applied current I app . Inward sodium current is not included here as it does not impact EAD generation due to the inactivation of this current at depolarized membrane potentials [41]. The calcium current activates instantaneously as a function of voltage, d ∞ (V ). Inactivation of the calcium current is governed by the gating variable f with steady-state inactivation f ∞ (V ) and time constant τ f . Activation of the potassium current occurs on a slower time scale and is described by the gating variable x with steady-state activation x ∞ (V ) and time constant τ x . The specific membrane capacitance is C = 1 µF/cm 2 .
The voltage-dependent activation and inactivation functions are sigmoids given by for y = d, f , and x, with half-(in)activation voltages θ y and slopes proportional to 1/σ y .
All of the parameters for the potassium current, except the maximal conductance g K , were kept the same as in the Sato model: τ x = 300 ms, reversal potential E K = −80 mV, and activation kinetics θ x = −40 mV, σ x = −5 mV. For the calcium current, we set τ f = 80 ms as in the Sato model, and then fit the remaining parameters to the voltage-clamp data from Chen et al. [9] shown in Figure 1A. In these recordings, Chen et al. measured the L-type calcium current in cardiomyocytes isolated from guinea pigs housed under 12h:12h light:dark cycles, with the lights turned on at Zeitgeber time 0 (ZT0, 7:00 AM) and turned off at ZT12 (7:00 PM). They performed voltage-clamp experiments in the morning (ZT 3) and at night (ZT 15), in which they held cardiomyocytes at -80 mV and then depolarized them in 10 mV increments from -70 to +60 mV. They found that at both times of day, the largest calcium currents were evoked at the +10 mV voltage step. Furthermore, the current density at this voltage step was significantly larger in the morning than at night.
We digitized their published I-V curves and normalized the data at ZT 3 and ZT 15 using the peak current density at +10 mV for each time of day. We then averaged these normalized curves to obtain a single curve to use as input to a parameter estimation algorithm. Specifically, we used an unconstrained nonlinear optimization routine (the Nelder-Mead algorithm fminsearch in MATLAB) and voltage-clamp simulations to find the parameter values E Ca = 60 mV, θ d = −7.3 mV, σ d = −8.6 mV, θ f = −13.3 mV, and σ f = 11.9 mV. These parameter values minimized the squared error between the model-generated I CaL I-V curve and the average normalized I-V curve from the voltage-clamp data. With the reversal potential and gating kinetics held fixed, we varied the maximal conductance and found that g CaL = 0.3 mS/cm 2 produced a model I-V curve (blue curve in Fig. 1B) similar to the experimental data from ZT 3, whereas the g CaL = 0.15 mS/cm 2 curve (red curve in Fig. 1B) was similar to the data from ZT 15.
To determine the maximal conductance g K , we performed current-clamp simulations and compared the model traces to the current-clamp recordings from Chen et al. [9] shown in Fig. 1C. For these simulations, we set I app = 0 and simulated for 10 ms with initial conditions V 0 = −80 mV, f 0 = f ∞ (−80) = 0.9963, and x 0 = x ∞ (−80) = 0.0003354. At t = 10 ms, we instantaneously set V = 0 to mimic the effect of a stimulating current pulse [34]. We then measure the action potential duration at 90% repolarization (APD90), which is the amount of time it takes for the voltage to return to 90% of its value before the spike. We find that setting g K = 0.1 mS/cm 2 with g CaL = 0.15 mS/cm 2 (blue curve in Fig. 1D) yields an APD90 for the model that is similar to the APD90 in the experimental data at ZT 15. Furthermore, setting g K = 0.1 mS/cm 2 with g CaL = 0.3 mS/cm 2 (red curve in Fig. 1D) gives a model APD90 that is very similar to the experimental APD90 at ZT 3.
Simulations of the single-cell model were performed using MATLAB R2017a (The Mathworks Inc., Natick, MA) and ode15s, a variable-step, variable-order solver for stiff ODEs.

Bifurcation analysis
The single-cell model was analyzed by decomposing it into a fast subsystem (Eq. 1-2) and a slow subsystem (Eq. 3), as in [66]. We then treat the slow variable x as a bifurcation parameter and study the bifurcation structure of the fast subsystem: The linearization of this system at a steady state (V * , f * ) is given by the Jacobian matrix At a steady state, we have that F (V, f ) = G(V, f ) = 0. To find steady states, we set f * = f ∞ (V * ) to satisfy G(V * , f * ) = 0, and then solve F (V * , f * ) = 0 for V * : Hopf bifurcation occurs when trace(J) = 0 and determinant(J) > 0. Saddle-node bifurcation occurs when trace(J) = 0 and determinant(J) = 0. Takens-Bogdanov (TB) bifurcation occurs when Hopf and saddle-node bifurcation points coalesce and the Jacobian matrix has two zero eigenvalues [14]. The conditions for this codimension-2 bifurcation are trace(J) = 0 and determinant(J) = 0, that is To find TB points, we simultaneously solve Eqs. 5-7. Bifurcations were also identified using the dynamical systems software package XPPAUT [16].

Spatial model
In cardiac tissue, neighboring cells are electrically coupled through gap junctions. The spatiotemporal evolution of the cellular membrane potential in a 2-D domain can be described by the following reaction-diffusion partial differential equation (PDE): where G x and G y are longitudinal and transverse conductances associated with the diffusion terms representing intercellular currents. To simulate a 2-D sheet of guinea pig cardiac tissue, we modified the monodomain reactiondiffusion MATLAB code developed by Hammer [26]. We solved the PDEs numerically on a 128 × 128 isotropic (G x = G y = 25 nS) grid using a finite-difference scheme for spatial derivatives, the explicit Euler method for time derivatives, and Neumann We used an S1-S2 cross-field stimulation protocol, with the first stimulus (S1) delivered to the left boundary of the domain at t = 0 with strength I app = 500 µA/cm 2 and a duration of 2 ms. The second stimulus (S2), was delivered to the bottom domain boundary at t = 810 ms with the same strength as S1 and a duration of 3 ms. This stimulation procedure generates spiral waves in our 2-D domain as shown in Fig. 6.

Elevated L-type calcium current in the morning can induce EADs
To investigate the role that circadian rhythmicity of the L-type calcium current plays in the electrical activity of guinea pig cardiomyocytes, we simulated an electrophysiological model of these cells (Eq. 1-3) with maximal conductance values corresponding to morning and evening time points: specifically g Ca = 0.3 mS/cm 2 at ZT 3 and g Ca = 0.15 mS/cm 2 at ZT 15. We determined these parameter values, along with the gating kinetics of the calcium current (Eq. 4), by fitting voltage-clamp data from Chen et al. [9] as described in the Methods section (see Fig. 1, A and B). This model, with maximal potassium conductance g K = 0.1 mS/cm 2 , can reproduce the circadian variation in action potential duration observed in current-clamp recordings (see Fig. 1, C and D). In the current-clamp data, APD90 is 11.5% greater at ZT 3 (228.0 ms) than ZT 15 (204.5 ms); in the model, APD90 is 16% greater at ZT 3 (225.5 ms) than ZT 15 (194.4 ms). The model enables us to explore the interaction between the potassium conductance and circadian variation of L-type calcium. We find that if g K is lowered to 0.05 mS/cm 2 , then the difference between morning and evening becomes more pronounced, with APD90 being 97.5% greater at ZT 3 (579.9 ms) than ZT 15 (293.6 ms), see Fig. 2. Moreover, the action potential at ZT 3 now exhibits secondary voltage depolarizations during the repolarization phase, known as early afterdepolarizations (EADs).

Dynamics of EAD generation
To understand the dynamical mechanism underlying the generation of these EADs, we follow Tran et al. [66] and perform a fast-slow decomposition of our model. As described in the Methods section, we study bifurcations in the fast (V, f ) subsystem (Eqs. 1-2) treating the slow variable x as a bifurcation parameter. The fast subsystem generally has three fixed points for small values of x and one fixed point for large values of x, forming a Z-shaped curve in the V, x plane (Fig. 3). The curve consists of an upper branch of depolarized fixed points, a middle branch of unstable fixed points, and a lower branch of hyperpolarized stable fixed points. As x is increased, the fixed points on the upper branch change from stable (solid red curve) to unstable (dashed black curve) at a subcritical Hopf bifurcation, where unstable limit cycles (open green circles) are born. These unstable periodic solutions are terminated at a homoclinic bifurcation with a saddle point on the middle branch of fixed points. As x is increased further, the upper and middle branches of fixed points approach each other and eventually coalesce, destroying these fixed points in a saddle-node bifurcation. With g CaL = 0.15 mS/cm 2 and g K = 0.1 mS/cm 2 , the Hopf and saddle-node bifurcations occur at passes through the V, x plane to the right of these values (x > x SN ), and therefore repolarizes monotonically without EADs. If g CaL is increased to 0.3 mS/cm 2 with g K held fixed (Fig. 3B), the Hopf and saddle-node bifurcation points move to the right (x HB = 0.376, x SN = 0.551), but so does the AP trajectory; here the AP repolarizes through the region x HB < x < x SN without EADs. Similarly, if g CaL is held fixed at 0.15 mS/cm 2 but g K is reduced to 0.05 mS/cm 2 (Fig. 3C), the trajectory repolarizes without EADs through the region x HB = 0.404 < x < x SN = 0.551. However, if g CaL is increased to 0.3 mS/cm 2 and g K is reduced to 0.05 mS/cm 2 , the model does exhibit EADs. The AP trajectory now repolarizes through a region where x < x HB = 0.753 and the fast subsystem contains a stable fixed point (Fig. 3D). This fixed point is a stable focus, and the trajectory exhibits damped oscillations until it reaches x HB . To the right of x HB , the fixed point is now an unstable focus and the trajectory exhibits one more voltage peak before repolarizing fully.

Takens-Bodganov bifurcation as an organizing center
Several different dynamical mechanisms can give rise to secondary oscillations that grow in amplitude, which is the EAD pattern typically observed in experiments. In the first mechanism to be characterized, stable limit cycles with growing amplitudes emerge from a supercritical Hopf bifurcation in the fast subsystem [66]. More recently, Kügler [34] demonstrated that EADs with growing amplitudes can also arise either from a delayed subcritical Hopf bifurcation in the fast subsystem, or along the unstable manifold of a saddle-focus fixed point in the full system. In the latter case, there is no Hopf bifurcation. The EADs explored earlier in this paper are of the subcritical Hopf type; see the bifurcation diagrams in Fig. 3. In these diagrams, the Hopf bifurcations occur relatively near a saddle-node bifurcation. This suggests that by varying another parameter in conjunction with the bifurcation parameter x, the Hopf and saddle-node bifurcations can be made to coalesce in a Takens-Bogdanov (TB) bifurcation. Indeed, Fig. 4A shows that the Hopf and saddle-node bifurcation points approach and collide with each other as x and g CaL are decreased simultaneously, with the TB bifurcation occurring at x = 0.411, g CaL = 0.0224 mS/cm 2 . We simulated the model with parameters chosen near the TB bifurcation point (g CaL = 0.02, g K = 0.01, θ x = −40, τ x = 1100) and observed EADs as shown in Fig. 4B-C. The eigenvalues of the full system linearized at the fixed point (V * , f * , x * ) = (−12.74, 0.4882, 0.3664), which corresponds to a location near the saddle-node bifurcation point shown in Fig. 4B, are λ 1,2 = 0.0025 ± 0.0066i and λ 3 = −0.0068. Thus, this fixed point is classified as a spiral saddle of index 2 [27]. The EADs arise due to the spiraling movement of the trajectory caused by the unstable manifold of the saddle focus, which is the second EAD-generating mechanism found by Kügler [34]. In this way, knowledge of the TB bifurcation can help identify parameter sets that produce different types of EAD dynamics.

Circadian variation of calcium and potassium currents
The voltage-clamp experiments of Chen et al. [9] revealed a day/night difference in L-type calcium current, with larger currents in the morning (ZT 3) than at night (ZT 15). Correspondingly, they observed longer duration action potentials in their current-clamp recordings at ZT 3 than at ZT 15. We used our model to simulate action potentials across the circadian cycle by assuming that g CaL follows a sinusoidal waveform with a period of 24 hours (Eq. 8), with a maximum of g CaL = 0.3 at ZT 3 and a minimum of g CaL = 0.15 at ZT 15. With g K held fixed at 0.1, we find that EADs do not occur at any time of day. Figures 5A and 4B show the APD90 values across a circadian cycle with g K = 0.1 and g K = 0.05 respectively. With g K = 0.05, we find that EADs occur over a large portion of the day (∼8 hours), specifically from ZT 23 to ZT 7. We simulated with two different g K values to represent the heterogeneity in potassium channel expression that has been found across different cells of the ventricular myocardium [69] or among different guinea pigs. However, it is plausible that g K may also be under circadian control, since it has been shown that IKs is modulated by K + channel interacting protein 2 (KChIP2) and that KChIP2 expression follows a circadian rhythm in mouse ventricles. In particular, co-expression of KChIP2 reduces IKs [39], and KChIP2 expression is higher at night than in the morning [29]. Thus we modeled a circadian rhythm in potassium channel conductance by assuming g K follows a sinusoidal waveform with a period of 24 hours (Eq. 9), with a maximum of g K = 0.1 at ZT 14 and a minimum of g K = 0.05 at ZT 2. Figure 5C shows the APD90 values across a circadian cycle with rhythms in both g CaL and g K . We find that EADs occur over an ∼5-hour portion of the day, from ZT 1 to ZT 6.

EADs lead to pathological electrical activity in 2-D tissue simulations
To explore whether the single-cell EADs triggered by circadian variation of ion channel conductances leads to abnormal electrical activity in cardiac tissue, we simulated a 2-D spatial domain as described in the Methods section. An S1-S2 stimulation protocol triggered spiral waves at both ZT 3 (g CaL = 0.3 mS/cm 2 ) and ZT 15 (g CaL = 0.15 mS/cm 2 ) with either normal (g K = 0.1 mS/cm 2 ) or low (g K = 0.05 mS/cm 2 ) potassium conductance (Fig. 6). Of these 4 scenarios, only the ZT 3 low g K combination exhibited EADs in the spatial model (Fig. 7). In addition, this combination led to the steepest APD restitution curve (Fig. 8), a commonly used indicator of the propensity for ventricular tachyarrhythmias [22,30,49]. To test this propensity, we next simulated heterogeneity in potassium channel conductance across the tissue with the leftmost 80% of the domain set to g K = 0.05 mS/cm 2 and the rightmost 20% set to g K = 0.1 mS/cm 2 . At ZT 15, the solution consists of a single spiral wave (Fig. 9A). However, at ZT 3, multiple spiral waves are born and collide into each other (Fig. 9B). This type of spiral wave break-up has been associated with ventricular fibrillation.

Cardiac arrhythmogenesis and circadian rhythms
Epidemiological studies have shown that the occurrence of life-threatening cardiac arrhythmias, such as ventricular tachycardia and ventricular fibrillation, exhibits time-of-day dependence with a peak in the morning hours [52]. For example, episodes of ventricular tachyarrhythmias recorded in patients with implantable cardioverter defibrillators were significantly increased between 8:00 and 11:00 AM [32]. Controlled laboratory studies indicate that the timeof-day fluctuations in adverse cardiovascular events are not solely due to daily rhythms in behavior and the external environment, suggesting that internally generated circadian oscillations influence cardiac arrhythmogenesis [8]. Normal electrical properties of the heart, such as electrocardiogram waveforms and heart rate, also demonstrate robust circadian rhythms [15,13]. The circadian system could exert this influence through two primary mechanisms: (1) regulation of cardiac electrophysiology by the central circadian clock in the hypothalamus through neurohumoral factors and the autonomic nervous system, or (2) local circadian clocks in cardiomyocytes themselves driving circadian rhythms in ion channel expression [3]. In this paper we considered the latter mechanism and explored how circadian rhythms in calcium and potassium conductances affect ventricular myocyte electrical activity across the day/night cycle.

Local cardiac circadian clock
Circadian clocks have been found in mammalian tissues throughout the body, including the heart. These peripheral clocks operate using the same molecular machinery as the central clock in the suprachiasmatic nucleus (SCN). The basic mechanism is a negative feedback loop in which the protein products of the clock genes Per and Cry inhibit their own production by repressing their transcriptional activator complex CLOCK-BMAL1. The timescales of the biochemical processes involved in this transcription-translation feedback loop lead to oscillations in the abundance of PER and CRY proteins with a period of approximately 24 hours. The expression of many other genes and proteins that are not necessarily integral to the clock mechanism are also influenced by CLOCK-BMAL1 and exhibit ∼24hour oscillations. Such clock-controlled genes (CCGs), including those that encode ion channels, can then modulate cellular processes in a time-of-day-dependent manner [76]. Oscillations in the expression of core circadian clock genes have been observed in the intact heart, cultured myocardial tissue, and isolated cardiomyocytes [3]. Studies in mice with cardiomyocyte-specific CLOCK mutations (CCM) and BMAL1 knockouts (CBK) demonstrate that 10% of the cardiac transcriptome is regulated by local circadian clocks in the heart [4,75]. Through these CCGs, the cardiomyocyte circadian clock impacts a variety of key cellular functions, including cardiac metabolism, signal transduction, contractility, and electrophysiology [45].

Circadian transcription of cardiac ion channels
Several ion channel subunits exhibit circadian rhythms in expression within the ventricles of animal models [3]. The levels of transcripts associated with Na + current (Scna5, Nav1.5, I N a ) [59], L-type Ca 2+ current (Cacna1c and Cacna1d, Cav1.2 and Cav1.3, I CaL ) [9,4], transient outward K + current (Kcnd2, Kv4.2, I to ) [65], ultra-rapidly activating delayed rectifier K + current (Kcna5, Kv1.5, I Kur ) [74], rapidly activating delayed rectifier K + current (Kcnh2, Kv11.1, I Kr ) [58], two-pore background K + current (Kcnk3, K2p3.1, I K2p ) [65], and gap junction current (Gja5 and Gja1, connexins Cx40 and Cx43) [64] oscillate over a 24-hour period. In some cases, rhythms in channel subunit gene expression have been related to day/night differences in electrophysiological properties and cardiac pacemaking. For example, elevated Kcna5 and Kcnd2 protein levels at ZT6 and ZT18, respectively, correlate with increased steady-state currents for I to and I Kur at those time points [74]. Potassium Channel Interacting Protein 2 (KChIP2), a regulator of I to , has been implicated in the circadian rhythm of cardiac repolarization. Jeyaraj et al. [29] showed that Kruppel-like factor 15 (Klf15) is a CCG that directly regulates KChIP2 expression, and that deletion of Klf15 abolishes the circadian rhythm in QT interval and increases suspectibility of mice to ventricular arrhythmias. However, Gottlieb et al. [23] concluded that KChIP2 is not a mechanistic link between the cardiac circadian clock and ventricular repolarization and arrhythmogenesis, based on their finding that KChIP2-deficient mice still have a circadian rhythm in QT interval. Rather, they suggest that Klf15 expression controls the transcription of other genes responsible for the circadian rhythm in repolarization and susceptibility to arrhythmias.

Circadian variation of L-type Ca 2+ current
In this paper we focused on circadian regulation of L-type Ca 2+ channels, due to the evidence supporting local cardiac clock control of these channels and the importance of L-type current for cardiac pacemaking. The α1D subunit of the L-type channel (Cacna1d) shows circadian variation in both mRNA and protein expression levels in the hearts of wild-type mice, but not in the hearts of CCM mice [4]. In guinea pigs, the α1C subunit of the L-type channel (Cacna1c) is rhythmic at the protein level with a peak at ZT3, which correlates to larger L-type calcium current at that time point [9]. Voltage-gated L-type Ca 2+ channels have also been proposed as a link between circadian oscillations in electrical activity and the molecular clock in SCN neurons [51,48,11,2] and retinal photoreceptors [31].
Although circadian variation of potassium channel expression has been observed in mouse and rat ventricles, the voltage-clamp studies of Chen et al. [9] did not find a significant time-of-day dependence for the major outward potassium currents (I Ks and I Kr ) in guinea-pig ventricular myocytes. Thus, in most of this paper we assume the potassium current to be constant throughout the day/night cycle. Instead, we consider the effect of circadian variation in  I CaL in the presence of lower or higher levels of I Ks , reflecting the heterogeneity in potassium channel expression one might expect to find across different cells or individuals. Nonetheless, in one set of simulations we do explore the effect of antiphase circadian variation of I CaL and I Ks (see Fig. 5).

Mathematical analysis of EADs
Mathematical modeling studies have shown that increased inward calcium current and decreased outward potassium current can elongate the cardiac AP and produce the pathological voltage oscillations known as early afterdpolarizations (EADs) [53,68,37,28]. To understand the genesis of EADs, minimal models of the cardiac AP have been analyzed using dynamical systems tools such as slow-fast decomposition and bifurcation theory. Tran et al. [66] showed that EADs involve supercritical Hopf and homoclinic bifurcations in the fast subsystem, and claimed that under periodic pacing the homoclinic bifurcation leads to chaotic behavior. Sato et al. [57] argued that deterministic chaos, rather than random fluctuations due to noise, is the primary cause of the irregular EAD dynamics frequently seen in cardiac experiments [57]. Kügler [34] showed that EADs can also arise from alternative dynamical mechanisms, such as delayed subcritical Hopf or limit point bifurcations of the fast subsystem. Furthermore, Kügler et al. [35] argued that a cascade of period doubling bifurcations underlies EAD chaos in both periodically paced and unpaced cardiomyocytes. These studies all decomposed the full model into fast and slow subsystems with a single gating variable in the slow subsystem. Kügler et al. [36] performed a slow-fast decomposition with two slow variables, and proposed a folded-node singularity of the slow flow as a novel mechanism for EAD generation. Vo and Bertram [70] analyzed the same model treating two variables as slow and also attributed EADs to folded-node singularities and their associated canard orbits. They demonstrated that the appearance of dynamical chaos under periodic pacing can be understood using the theory of candard-induced mixed-mode oscillations [5].
In this paper, we utilized the same 3-variable model for cardiac APs introduced in [57] and studied in [34,70], but we refit the parameters of the L-type calcium current to the voltage-clamp data of [9]. With these parameters, when the model is analyzed with a 1-slow-2-fast structure, the EADs are generated by a subcritical Hopf bifurcation in the fast subsystem. This is one of the EAD mechanisms explored in [34]. We then showed that a Takens-Bogdanov (TB) bifurcation is present in this model, and that near the TB point we can find EADs generated by the unstable manifold of a saddle-focus fixed point of the full system. This is the other EAD mechanism explored in [34]. Thus, the TB bifurcation serves as an organizing center for the dynamics and helps connect some of the different types of EADs that have been observed previously.

Modeling of cardiac tissue
To study how circadian variation of ionic conductances affects cardiac excitability at the tissue level, we simulated a 2-D spatial model using reaction-diffusion PDEs and an S1-S2 stimulation protocol. The spatial model exhibited spiral wave solutions at both circadian time points (ZT3 and ZT15) and with both low and high potassium conductance (g K = 0.05 and 0.1). Under the conditions where the single-cell model exhibits EADs (ZT3 with g K = 0.05), the spiral waves in the spatial model had a faster propagation speed, analogous to the heart beating faster as during ventricular tachycardia. When spatial heterogeneity in potassium conductance was introduced, the time of day where the single-cell model exhibits EADs produced spiral wave break-up, a behavior associated with ventricular fibrillation [17]. It is generally accepted that EADs at the cellular level can lead to arrhythmias, such as polymorphic ventricular tachycardias (PVT) and Torsade de Pointes (TdP), at the tissue level [70]. Modeling studies have shown that single-cell EADs can cause wave initiation, that these EADs can synchronize to form 2D wave patterns, and that meandering waves in heterogeneous tissue can give rise to the classic ECG appearances of PVT and TdP [71,12,7]. Vandersickel et al. [68] performed a systematic study of single-cell EAD excitations and their 2D manifestations in a model of human ventricular tissue. However, there are still many open questions about how EADs progress to perpetuating arrhythmias [67].

Conclusion and future directions
The main finding of this paper is that circadian rhythms in L-type calcium conductance can lead to early afterdepolarizations at certain times of the day in a model of guinea pig ventricular myocytes. To our knowledge, this is the first study to consider how the cardiomyocyte circadian clock influences the genesis of EADs. We propose that circadian rhythms in EAD occurrence may contribute to the time-of-day-dependent patterns observed in ventricular tachyarrhythmias and sudden cardiac death. However, to establish this connection there are some limitations of our study that would need to be addressed, as discussed below.
First, Zeitgeber times (ZT) for the guinea pig experiments need to be related to real-world time for humans. Guinea pigs are a commonly used animal model for cardiac electrophysiology, as the shape of guinea pig action potentials are more similar to human APs than are the APs of smaller rodents such as mice. On the other hand, guinea pigs are not a commonly used animal model for circadian experiments, and they don't have particularly strong sleep/wake rhythms [10]. Guinea pigs are crepuscular, meaning they are most active at dawn and dusk and are neither nocturnal nor diurnal [38]. To extrapolate results from studies with nocturnal animal species to diurnal humans, a phase shift of approximately 12 hours is often assumed. For example, in a study assessing the relevance of circadian rhythms for administration of chemotherapy, 7 hours after light onset (ZT7) for mice was taken to be "middle of the night" for humans [24]. However, a recent study found that many cardiovascular drug targets exhibit circadian rhythms in gene expression with a similar phase relationship in mouse and human heart tissue, including several genes encoding L-type Ca 2+ channel subunits [55]. In our guinea pig simulations, EADs occurred between ZT0 and ZT5 (Fig. 5C). Assuming a similar phase relationship between guinea pigs and humans, this corresponds to an increased likelihood of EAD-induced arrhythmias in the first few hours after waking up, in accordance with the peak time window for sudden cardiac death found in epidemiological studies [47,72,63].
Second, in this study we employed a minimal model of cardiac AP generation consisting of a single Ca 2+ current and a single K + current. The advantage of this approach is that the low dimensionality of the model facilitates bifurcation analysis and an understanding of how circadian rhythms affect the dynamics underlying EAD generation. A disadvantage is that the model is lacking descriptions of some specific types of ion channels that may be relevant for daily variation of cardiac electrical properties. For example, Nav1.5 sodium channels and Kv11.1 (mERG) potassium channels display circadian rhythms in transcription in mouse hearts [59,58]. Moreover, cardiomyocyte-specific deletion of Bmal1 abolishes circadian oscillations in both of these channels, suggesting they are under the control of the local cardiac molecular clock [3]. In future work, we plan to investigate how circadian rhythms in these conductances affect the propensity for EADs and arrhythmias in more detailed models of cardiac electrophysiology that include many types of ionic currents and intracellular calcium dynamics.
Due to the critical role that I CaL plays in EAD formation, L-type Ca 2+ channels have been identified as a promising therapeutic target for suppressing EADs and their arrhythmogenic consequences [43,42,44]. Based on the results of our study, we suggest that special attention should be paid to the time of day that drugs targeting L-type channels for EAD suppression are taken in order to enhance their effectiveness. Tailoring the timing of drug administration based on circadian factors, known as chronomedicine or chronopharmacology, is an emerging area of precision medicine with clinical trials showing dosing-time-dependent efficacy or toxicity across several conditions, including hypertension and other cardiovascular disorders [54,6]. Computational modeling of how the circadian clock modulates therapeutic targets can be used to help predict the optimal dosage time to maximize efficacy or minimize side effects [1].
Cardiotoxicity is the leading cause of drug development discontinuation and withdrawal of drugs from the market [54,18]. Multiple drugs that have been pulled from the market for causing fatal TdP have the unintended side effect of blocking Kv11.1 (hERG) potassium channels, and screening for ERG block is now a mandatory requirement for new pharmaceuticals [46]. ERG block is a sensitive but not specific measure of TdP risk, i.e. it gives few false negatives but false positives may be preventing safe drugs from entering the market [50]. The Comprehensive in vitro Pro-arrhythmia Assay (CiPA) is a new global initiative to create guidelines for the assessment of drug-induced TdP that recommends a central role for computational modeling of ion channels and in silico evaluation of compounds [19,40]. As noted above, many cardiac ion channels exhibit circadian oscillations, including ERG. Thus, we propose that circadian clock modeling should be incorporated into the CiPA paradigm for assessing drug-induced cardiotoxicity.