A Mathematical Model of Neonatal Rat Atrial Monolayers with Constitutively Active Acetylcholine-Mediated K+ Current

Atrial fibrillation (AF) is the most frequent form of arrhythmia occurring in the industrialized world. Because of its complex nature, each identified form of AF requires specialized treatment. Thus, an in-depth understanding of the bases of these arrhythmias is essential for therapeutic development. A variety of experimental studies aimed at understanding the mechanisms of AF are performed using primary cultures of neonatal rat atrial cardiomyocytes (NRAMs). Previously, we have shown that the distinct advantage of NRAM cultures is that they allow standardized, systematic, robust re-entry induction in the presence of a constitutively-active acetylcholine-mediated K+ current (IKACh-c). Experimental studies dedicated to mechanistic explorations of AF, using these cultures, often use computer models for detailed electrophysiological investigations. However, currently, no mathematical model for NRAMs is available. Therefore, in the present study we propose the first model for the action potential (AP) of a NRAM with constitutively-active acetylcholine-mediated K+ current (IKACh-c). The descriptions of the ionic currents were based on patch-clamp data obtained from neonatal rats. Our monolayer model closely mimics the action potential duration (APD) restitution and conduction velocity (CV) restitution curves presented in our previous in vitro studies. In addition, the model reproduces the experimentally observed dynamics of spiral wave rotation, in the absence and in the presence of drug interventions, and in the presence of localized myofibroblast heterogeneities.


Introduction
Atrial fibrillation (AF) is the most frequent form of arrhythmia occurring in patients with or without the risk for heart failure [1][2][3][4][5][6][7][8]. Each identified form of AF (paroxysmal, persistent, or chronic) requires specialized treatment. Therefore, an in-depth understanding of the underlying mechanisms is required for optimal management in the treatment of such arrhythmias.
Chronic AF is known to induce functional tissue remodeling through changes in ion channel expression and activity in the cardiomyocyte membranes. Of the various ion channels that have been identified in atrial cardiomyocytes, one family of particular interest is the Kir3.x. These channels produce an inwardly-rectifying K + current (I KACh ) that is normally controlled by vagal nerve activity through its transmitter acetylcholine (ACh). There is growing consensus that during chronic AF, I KACh becomes constitutively active (I KACh-c ) [9][10][11] and hyperpolarizes the cell membrane. This leads to abbreviation of the action potential duration (APD), and stabilization of reentrant electrical circuits that characterize AF [1]. However, despite decades of research on a variety of mammalian hearts, the role of ion-channel remodeling in AF, remains incompletely understood.
Development of a better biophysical understanding of ion channel remodeling in AF requires thorough in vivo investigations in the human heart. However, a more practical approach involves studying AF in cell cultures in vitro, prior to ex vivo and in vivo trials. Primary cultures of neonatal rat atrial cardiomyocytes (NRAMs) have been used by researchers interested in cell signaling [11][12][13]. Recently, we have demonstrated that NRAMs can generate I KACh-c , in primary cultures, thereby making them attractive sources for arrhythmogenic substrates that can be used to study chronic AF [11].
Although in vitro and in vivo studies have often provided significant insights into the mechanisms underlying reentry initiation in both healthy and diseased cardiac tissue, they have a few limitations. One cannot easily address complex questions related to electrophysiology, that would involve challenging techniques like patch-clamp in parallel with optical mapping measurements, or sub-surface measurements of electrical signals in intact cardiac tissues. Furthermore, independent experimental interventions involving gradual modifications of cell properties, such as gap junctional conductance, coupling between different cells types, localized changes in conductivities of specific ionic channels etc. are very difficult to obtain in an experiment. This is where computer models find wide-scale applicability, from basic science [14][15][16][17][18][19][20][21][22][23] to direct clinical applications [24][25][26][27][28][29][30][31][32][33].
The basis for any computational study in cardiac tissue is a model for the cardiac cell. However, a model for the NRAM is still lacking. Development of this model is important, because it would allow researchers to perform specific in silico studies next to in vitro experiments, thereby increasing the possibilities for the development of biological insight into arrhythmias like AF.
In this study, we formulate the first complex ionic model for NRAM monolayers, based on patch-clamp data from recent literature. To remain consistent with in vitro experiments, we incorporate natural cellular heterogeneity, and a baseline 15-20% randomly distributed myofibroblasts into these monolayers. We then generate spiral waves in these monolayers by burst pacing (BP) and reproduce spiral wave dynamics as seen in experiments, both under normal conditions and under pharmacological inhibition. In particular, we simulate the effect of the atrial-specific drug tertiapin-Q, a blocker of I KACh , on BP-induced arrhythmias, as reported by Bingen et al. [11] In addition, we reproduce spiral wave dynamics in a monolayer with localized myofibroblast heterogeneity.

Formulation of the single cell ionic model
General. The mathematical model of the native NRAMs was developed from an electrical circuit of the cell membrane. It consists of a membrane capacitance, in parallel with voltagedependent conductance branches that include batteries, and with pumps and exchangers. Thus, the electrophysiological behavior of a single cardiomyocyte can be described using the following ordinary differential equation, where V is the transmembrane voltage in millivolts (mV), t is time in milliseconds (ms), I ion is the total transmembrane ionic current density in microamperes per square centimeter (μA/ cm 2 ), I stim is the externally applied stimulus current, and C m is the cell capacitance per unit surface area in microfarad per square centimeter (μF/cm 2 ). In our simulations with the single cell, we applied a stimulus current of strength 7 pA for a duration of 5 ms. I Na represents the fast Na + current, I CaL , the L-type Ca 2+ current, I K1 , the inward-rectifier K + current, I to , the transient outward K + current, I CaT , the T-type Ca 2+ current, I Cab , the background Ca 2+ current, I NCX , the Na + /Ca 2+ exchanger current, I NaK , the Na + /K + pump current, I f , the hyperpolarization-activated funny current, I Nab , the background Na + current, and I KACh , the acetylcholine-mediated K + current. All conductances G X were measured in nanosiemens per picofarad (nS/pF), while intracellular and extracellular ionic concentrations ([X] i , [X] o ) were expressed in millimoles per liter (mM). The detailed equations for each of the ionic currents, are specified in the S1 Appendix of the supplementary material, along with a list of model parameters and initial values. The schematic diagram of the cell is shown in Fig 1. Membrane currents. Fast Na + current (I Na ): I Na was modeled according to the three-gate formulation first used by Beeler and Reuter [34].
Here m represents the activation gate, h, the fast inactivation gate, and j, the slow inactivation gate. The formulation for the channel kinetics was based on experimental data from Ramos-Mondragόn et al., [35] and Voitychuk et al. [36] We used the normalized activation (control) data from Fig 4D of Ramos-Mondragόn et al., [35] to construct the steady-state activation curve (m 1 ) of our model. For steady-state inactivation, we used the (control) data from Fig 5D of Voitychuk et al. [36] Fig 2A shows the steady-state activation and inactivation curves generated using our model, compared with the activation and inactivation kinetics from Ramos-Mondragόn et al. [35] and Voitychuk et al., [36] respectively. The time constants for the activation and slow inactivation gates are taken from the Luo-Rudy guinea-pig ventricular cell model, [37] and scaled appropriately for room temperature [38]. We retain these values for the time constants until more experimental evidence is available (Fig 2B and 2C). The time constant for fast inactivation was modeled on the basis of data derived from Fig 5E of Voitychuk et al., (Fig 2C) [36]. The maximal conductance (G Na ) was determined by fitting the maximum upstroke velocity (dV/dt max ) that we measured in vitro. This value was found to be in consonance with the recordings of Voitychuk et al. [36] Finally, current-voltage relationships (I-V curves) were constructed from the model formulation and fitted with experimental data from Fig 4B of Ramos-Mondragόn et al., (Fig 2D) [35]. The individual equations for these three gates governing Na + -channel kinetics, together with the time constants and initial values, are provided in the S1 Appendix of the supplementary material.
L-type Ca 2+ current (I CaL ): Although I CaL plays an important role in the generation of an AP in ventricular cells, in atrial cells its magnitude is comparatively less (cf. Figs 2B and 3B of Avila et al., [39]). Furthermore, the electrophysiological influence of L-type Ca 2+ in NRAMs is not well studied. Therefore, we used a more general formulation of I CaL , namely, the one according to the model by ten Tusscher et al. [40] It consists of a voltage-dependent activation gate d, a voltage-dependent inactivation gate f, and an intracellular calcium-dependent inactivation gate f Ca . The driving force was modeled Schematic diagram of the model of a cultured neonatal rat atrial cardiomyocyte (NRAM). The model includes fast Na + current (I Na ), background Na + current (I Nab ), L-type Ca 2+ current (I CaL ), T-type Ca 2+ current (I CaT ), transient outward K + current (I to ), background Ca 2+ current (I Cab ), hyperpolarization-activated funny current (I f ), inward-rectifier K + current (I K1 ), rapid and slow delayed rectifier K + currents (I Kr and I Ks ), ultra-rapid outward K + current (I Kur ), Na + /Ca 2+ exchanger (I NCX ), Na + /K + -ATPase (I NaK ), and a constitutivelyactive acetylcholine-mediated K + current (I KACh-c ). The sarcoplasmic reticulum (SR) is divided into 2 compartments for uptake (NSR) and release (JSR). Intracellular Ca 2+ is transported from the NSR to the JSR via diffusion (I tr ). In JSR, the Ca 2+ is buffered by calsequestrin (CSQN), and released into the cytosol (I rel ) via ryanodine receptors (RyR). While most of the intracellular Ca 2+ is pumped into the NSR via the SERCA (I up ), some fluxes passively diffuse back to the cytosol (I leak ). The Ca 2+ in the cytosol is buffered by troponin (TRPN) and calmodulin (CMDN). with a Goldmann-Hodgkin-Katz equation, ignoring the permeability for the Na + and K + ions. The steady-state activation and inactivation curves for the voltage-dependent gates (i.e., d 1 and f 1 ) were obtained by fitting the patch-clamp (control) data from the experiments of Avila et al, on NRAMs ( Fig 3A) [39]. Unfortunately, no literature was available for the activation and inactivation times of I CaL in NRAMs. Thus in our model, we retained the activation and inactivation time constants from the neonatal rat ventricular cardiomyocyte (NRVM) model by Korhonen et al. [41] (Fig 3B and 3C), with a 10 ms correction to the activation time constant (τ d ).
Korhonen et al. [41] obtained the value for Ca 2+ -independent inactivation from the experimental records of Pignier et al., [42] who used NRVMs for their measurements. The formulation for the [Ca] i -dependent inactivation was also taken from the NRVM model by Korhonen et al., [41] (Fig 3D). The maximal conductance for the L-type Ca 2+ channel was adjusted along with maximal conductances of the I K1 , till a reasonable fit of the control APD 80 restitution curve was obtained. The normalized I-V curve for I CaL , as obtained from the present model is shown in Fig 3G and overlaid on experimental data derived from Avila et al. [39] T-type Ca 2+ current (I CaT ): The T-type Ca 2+ current was formulated on the basis of Dokos et al., [43] as: where b represents the voltage-dependent activation gate and g, the voltage-dependent inactivation gate. The steady-state activation and inactivation curves were obtained by fitting the (control) data from the experiments of Avila et al., (Fig 3E) [39]. As in the case of I CaL , data on activation and inactivation times of T-type Ca 2+ channels are still lacking. Therefore, we used time constants (τ b and τ g ) from the NRVM model by Korhonen et al. (Fig 3F) [41] The maximal conductance (G CaT ) was adopted from the NRAM data of Avila et al. [39] The normalized I-V curve for I CaT generated using our model is shown in Fig 3H, together with I-V curve data derived from control experiments of Avila et al., on NRAMs [39]. Background currents (I Cab , and I Nab ), and hyperpolarization-activated funny current (I f ): Two background currents were directly taken from the adult rat ventricular myocyte model by Pandit et al. [44] These are linear (Ohmic) currents, representing the leakage of ions through non-specific ion channels in the NRAM-membrane. The maximal conductances of these are chosen such as to obtain the correct resting membrane potential, and are listed in the S1 Appendix of the supplementary material. There is substantial evidence regarding the presence of a hyperpolarization-activated, funny current (I f ) in NRVMs [45][46][47][48]. However, whether this current is also present in NRAMs is unclear. In line with the concept that it is present in both sino-atrial nodal cells and ventricular cells, we assume that it is also present in the atrial cells, but perhaps, not so significantly. Therefore, we include I f in the NRAM model, based on the formulation from the ventricular model by Korhonen et al. [41] The activation kinetic of the voltage-dependent channel y is shown in Fig 4A, together with the voltage-dependency of its time constant τ y (Fig 4B). We adopted the same value of the maximal ionic conductance for I f as Korhonen et al. [41] However, the presence or absence of this current did not make any significant difference to the morphology and characteristics of the AP or to the resting membrane potential.
Transient outward K + current (I to ): The transient outward K + current was modeled as: It comprises a voltage-dependent activation gate r, a voltage-dependent fast inactivation gate s, and a voltage-dependent slow inactivation gate s slow . A comparison of I-V curves for outward K + currents recorded from NRAMs and NRVMs show that the trends in I to are the same in both cell-types (cf. Fig 4B of Saygili et al., [49] and Fig 1B of Cahill et al. [48]), although the magnitudes are different (ventricular cells produce more I to than atrial cells). Thus, we used the same formalism for the steady-state activation (r 1 ) and inactivation curves (s 1 and s slow1 ), and the activation (τ r ) and inactivation time constants (τ s , and τ slow ), as in the rat ventricular cardiomyocyte model by Pandit et al. [44] The ion channel kinetics for the transient outward K + current are shown in Fig 5A-5C. Fig 5D shows the normalized I-V curve for simulated I to , overlaid on experimental data derived from Ramos-Mondragόn et al. [35] Sustained outward K + current (I Ksus ): In order to compare the I-V curve generated by the model, with I-V curve data derived from the control experiments of Ramos-Mondragόn et al., on NRAMs [35], the sustained outward K + current was modeled as a sum of three currents: the ultra-rapid current I Kur , the rapid-delayed rectifier current I Kr , and the slow-delayed rectifier The formulation for the first component I Kur was taken from the Bondarenko model [50]. The second and third components of I Ksus , namely I Kr and I Ks , were taken from the NRVM model by Hou et al. [51] All three currents are described in the S1 Appendix of the supplementary material. The maximal conductances for the three components (G Kur , G Kr , and G ks ) were adopted directly from the ventricular model by Hou et al. [51] However, the net sustained current I Ksus was adjusted (by the factor 0.16 in Eq 7) to reproduce the measured control APD 80 restitution curve. Steady-state activation and inactivation kinetics of the ionic conductances involved in producing the three component currents are illustrated in Fig 6A-6F. Fig 6G shows the normalized I-V characteristic curve for simulated I Ksus , overlaid on experimental data derived from Ramos-Mondragόn et al. [35] Inward-rectifier current and K + leak ( I K1 ): The NRAMs showed less inward-rectifier current than the NRVMs. For I K1 , the basic formulation of the model by Hou et al. [51] was used, with additional adjustments based on our patch-clamp recordings: The current consists of a [K] o -dependent part, a voltage-dependent inward rectification factor, and a linear leakage (last term). Our patch-clamp recordings of the inward-rectifier current were contaminated with a K + leakage current, which shifted the reversal potential of the net measured current by~10 mV. This shift was also observed in the records of Saygili et al. who Steady-state characteristics and time constants of the sustained K + current (I Ksus ) generated using the model. The sustained current was a composition of the ultra-rapid K + current (I Kur ), the rapiddelayed rectifier K + current (I Kr ), and the slow-delayed rectifier K + current (I Ks ). A, steady-state activation (ua 1 ) and inactivation (ui 1 ) curves of I Kur . B, activation (τ ua ) and inactivation (τ ui ) time constants of I Kur . C, steady-state activation (Xr 1 )and inactivation (Rr 1 ) curves of I Kr . D, activation time constant (τ Xr ) of I Kr . E, steady-state activation curves (Xs1 1 and Xs2 1 ) of I Ks . F, fast (τ Xs1 ) and slow (τ Xs2 ) inactivation time constants of I Ks . G, normalized I-V curve for sustained K + current (I Kur + I Kr + I Ks ), overlaid on patch-clamp data (open circles) derived from Ramos-Mondragόn et al. [35] doi:10.1371/journal.pcbi.1004946.g006 showed the complete I-V curve for the inward rectifier (I K1 ) in the NRAMs [49]. Thus, based on our patch-clamp data, the maximal conductance was adjusted to 47.5% of that of a NRVM [51], and the reversal potential for I K1 was shifted by +10 mV, ensuring that the total current I K1 was not a pure inward-rectifier, but a composite of the pure I K1 and the K + leak. An example of the inward rectifier current is shown in Fig 7B. Acetylcholine-mediated K + current (I KACh ): Unlike the situation in NRVMs, in NRAMs the acetylcholine (ACh)-mediated K + current (I KACh ) contributes to the repolarization reserve. We performed whole-cell voltage-clamp experiments, specific for the measurement of this current. We measured I KACh as the instantaneous time-dependent increase in current upon hyperpolarization (I KH ) [52][53]. Hyperpolarizing voltage clamp steps (5s) from -40 mV (to inactivate Na + current and T-type Ca 2+ current) evoked inward currents. As expected, the current exhibited pronounced inward rectification, with tiny current at potential positive to -50 mV ( Fig  7A). The addition of tertiapin-Q (100 nM) to the modified extracellular solution attenuated I KH (Fig 7B). Using our patch-clamp recordings of the time-dependent hyperpolarization-activated current I KH in Fig 7A and the inward-rectifier current I K1 in Fig 7B, we obtained the tertiapin-sensitive I KACh  The time-dependent hyperpolarization-activated current (I KH ) and its block by tertiapin. A, The voltage-clamp recordings depict large inward currents evoked by hyperpolarizations (5s), whereas, very little current was activated during steps to -50 mV and more positive potentials. In addition to pharmacological blocking, while performing voltage-clamp recording, Na + current (I Na ) and T-type Ca 2+ current (I CaT ) were inactivated by holding the cells at the potential of -40 mV. B, In the presence of tertiapin (100 nM), which selectively suppresses the constitutively active acetylcholine (ACh)-mediated K + current I KACh-c component of I KH , only inward rectifier current (I K1 ) is left with K + leak. C, Tertiapin-sensitive currents or I KACh-c (B subtracted from A). D, comparison of I-V characteristics, generated by our model (dashed line for control and solid line for tertiapin-Q treatment) and mean±SEM (open circles for control and filled circles for tertiapin-Q treatment) current density based on 8 neonatal rat atrial cardiomyocytes (NRAMs). E, representative AP-recordings from a single NRAM in control condition (black), upon treatment with 2 μM Carbachol (red), and upon subsequent treatment with 100 nM tertiapin-Q (gray). F, representative AP optical mapping traces recorded from a NRAM monolayer, under control condition (black), upon treatment with 2 μM Carbachol (red) and upon treatment with 100 nM tertiapin-Q (gray). G, compensated I-V characteristics for monolayer model, overlaid on the experimental single-cell I-V curves of D.
Near-overlap of the I-V curves in the physiological range of voltages indicates that in single cells I KACh is present, but not prominently active (Fig 7D). To verify this finding we performed additional current-clamp experiments. In particular, we evoked APs from 4 different NRAMs, in 4 independent control experiments, followed by subsequent treatment with 2 μM carbachol (the nonselective muscarinic receptor agonist). With this concentration of carbachol, I KACh is expected to be fully activated, thereby leading to abbreviation of APs in NRAMs. AP recordings, during carbachol-treatment, showed an average~35% decrease in APD 50 and an approxi-mate~40% decrease in APD 80 , confirming that I KACh was being abundantly activated in NRAMs by carbachol (Fig 7E, red line). Further treatment of these cell with 100 nM tertiapin-Q caused the original AP characteristics to be restored and lengthened of APD 80 (Fig 7E, gray  line). This finding suggests that isolated NRAMs expressed potentially active Kir3.x channels, which contribute to I KACh . However, this current is dormant on single NRAMs to produce substantial amounts of I KACh-c .
In monolayers, however, I KACh is found to be constitutively active and has been established to play a key role in the development of a proarrhythmic tissue substrate, because block of I KACh drastically increases APD, whereas further opening with 2 μM carbachol does not have any substantial effect (Fig 7F) [9][10]55]. Previously, we showed that complete block of the Kir3.x channels could lead to drastic changes in the overall electrophysiological behavior of the NRAM monolayer, such as near-tripling of APD 80 , reduction in APD 80 dispersion, and consequent termination of burst pacing-induced arrhythmias [11]. Therefore, we modified our formulation for the I KACh in higher dimensional models to make it constitutively active (I KACh-c ; Fig 7G). However, measurement of I KACh at physiological voltages in the presence of Carbachol requires an extracellular environment that contains abnormally high concentrations of K + , which would cause us to study a condition that will not reflect the actual situation in monolayers, which we aim to model and study in silico. Therefore we modified the single-cell formulation of this current, based on the concept derived from Dobrev et al (see Fig 1A of Dobrev et al.) [9], and increased the gap between the I-V curves for I K1 and I KACh in the physiological range of voltages. The modified formulation is shown in Eq 10.
Na + /Ca 2+ exchanger (I NCX ), and Na + /K + pump (I NaK ) currents: The Na + /Ca 2+ exchanger current was taken from the adult rat ventricular cardiomyocyte model by Pandit et al., [44] whereas the formulation for the Na + /K + pump current was adopted from the Luo-Rudy model of guinea pig ventricular cardiomyocyte [37], keeping the same Na + and K + concentration levels as used in our patch-clamp experiments (see Methods section). The formulations for these currents and their maximal conductances are listed in the S1 Appendix of the supplementary material.
Intracellular ion dynamics. For simplicity, the intracellular Na + and K + concentrations were kept constant in time. Dynamics of Ca 2+ changes was adapted from the NRVM model by Korhonen et al. [41] The sarcoplasmic reticulum (SR) was modeled as a composition of two compartments, the SR uptake or the network sarcoplasmic reticulum (NSR), and the SR release or the junctional sarcoplasmic reticulum (JSR) (Fig 1). The volume of JSR was assumed to be 10% of the total volume of the SR. [41] The Ca 2+ flux between the NSR and the cytosol through the SERCA, was adopted from the rabbit ventricular cardiomyocyte model by Shannon et al., [56] and modified in accordance with the NRVM model by Korhonen et al. [41] With the same assumption as Korhonen et al., [41] that the calsequestrin buffer was located close to the ryanodine receptors (RyRs), the complex intermolecular coupling within an RyR cluster was replaced by a single average RyR, and the opening of the RyR was modeled as being regulated by the cytosolic Ca 2+ concentration near the SR and the JSR. Our model accounted for Ca 2+ buffering through troponin and calmodulin.
The origins of the formulations of all the currents and the data used to adjust the model are summarized in Table 1.

Features of the single-cell ionic model
Single cell. Fig 8A inset shows an action potential (AP) generated using our model (dashed line), overlaid on an AP measured with whole-cell patch-clamp (solid line), at a membrane potential of~-85mV, which was established by a holding current to obtain a standard condition for the measurement of excitability. The AP shows the characteristic sharp upstroke, followed by a decay (repolarization phase), with complete absence of a notch and a dome. In the absence of a holding current, the model generated a resting membrane potential (RMP) of -72.0 mV, which was in consonance with the value -73.0 ± 2.5 mV, derived from the patchclamp measurements of Voitychuk et al. [36] dV/dt max was estimated to be 114 mV/ms, in line with 114.80±24.77 mV/ms measured in our experiments, and those of Voitychuk et al. [36]  APD 80 and APD 50 also showed good agreement with our experimental values (see Table 2). The calcium transients generated by the model are shown in Fig 9A and 9B. The mean concentrations of Ca 2+ in the uptake and release compartments were~775 μM and~715 μM, respectively, as opposed to their corresponding counterparts of~750 μM and~650 μM in ventricular cells [35], The mean free cytoplasmic Ca 2+ was found to be~0.36 μM in simulated atrial cells, as opposed to 0.475 μM in ventricular cells (Fig 8A of Korhonen et al. [41]). In general the amplitude of Ca 2+ cycles was found to be smaller than in neonatal rat ventricular cells [41]. Fig 10A-10F show the major ionic currents underlying the AP generated using our singlecell model. The fast Na + current (Fig 10A) generated by our model shows good agreement with the records of Ramos-Mondragόn et al., [35] (cf. Fig 4A of Ramos-Mondragόn et al., [35]) Unlike the ventricular model [41], our NRAM model shows reduced I CaL , but comparable I CaT .   ( Fig 10D) The calculated I to (Fig 10C) is also much smaller than in the ventricular model [41]. This could explain the absence of the notch in the morphology of the AP. The impure inward rectifier ( I K1 ) current in Fig 9C is a composite of the pure I K1 and a K + leakage current. Overall our model shows strong correspondence of the dynamical behavior of ionic currents with the available direct experimental measurements, and reproduces good integral cell properties such as action potential shape. We provide the numerical code for the single cell model in S1 Code. The next step is to verify wave propagation in two-dimensional (2D) preparations.

Extension of the single-cell ionic model to higher dimensions (2D monolayers)
General formulation. The monodomain formulation was used to extend our single-cell model to higher dimensions. Thus a 2D sheet of cardiac cells, representing a confluent monolayer, was modeled with the following partial differential equation: The second term on the right hand side of equation (Eq 11) takes care of the intercellular coupling. The coupling coefficient D is a symmetric, rank-3 tensor, whose elements determine the anisotropy of cardiac tissue. However, ordinary confluent in vitro monolayers are isotropic. Therefore, the 2D in silico studies presented here, made use of isotropic simulation domains where the diffusion tensor was replaced by a scalar diffusion constant. The temporal part of Eq 11 was integrated using the forward Euler method with time step δt = 0.02 ms, whereas, the spatial part was integrated using a centered finite-differencing scheme with a space step δx = δy = 0.00625 cm and Neumann (zero flux) boundary conditions.
Modeling cardiac monolayers. Our in vitro monolayers were cultured on circular coverslips in standard 24-well plates which have a basal diameter of 15.6 mm. These cultures showed an average 15-20% randomly-distributed myofibroblasts [11,[57][58][59][60]. Therefore, to simulate cardiac monolayers that closely resembled our in vitro cultures, a circular simulation domain of physical diameter 15.6 mm was used. The myofibroblasts were modeled according to the formulation of MacCannell et al., [61] which treats them as passive circuit elements that are electrically coupled to cardiomyocytes. Simulated cardiomyocytes, together with 15-20% simulated MacCannell-type myofibroblasts [61] were embedded on a rectangular grid, as described previously [22]. Natural intercellular variability was incorporated within the monolayer as described in Methods. The cells on the grid were then coupled diffusively (Eq 11). Fig 8A and 8B show a comparison of APs measured at 1Hz electrical pacing using optical mapping, with corresponding voltage traces generated from in silico monolayers. The optically mapped APs have similar APDs but aberrant rise and decay times because of filtering associated with the optical mapping technique.
Wave propagation in 2D. In monolayer control experiments, NRAMs exhibited a planar conduction velocity (CV) of~20-25 cm/s. We found that, to obtain such velocity, a diffusion coefficient of 0.00012 cm 2 /ms was required, which resulted in a CV of 22.2 cm/s. Optimal time and space steps for our model were δx = 0.00625 cm, δt = 0.02 ms, which resulted in an error of about 4% for the velocity of wave propagation and was well within tolerable limits [21]. Fig  11A shows baseline electrical propagation in an in vitro cardiac monolayer at 1Hz electrical pacing. The corresponding in silico activation map generated using our monolayer model is shown in Fig 11B. The two figures show excellent correspondence with roughly the same number of activation contours and similar isochronal line spacing.
In order to develop the model for I KACh-c in 2D, we fitted the APD 80 restitution (APDR) curves with experimental data from existing literature. In experiments, both APDR and CVR curves were measured via a dynamic protocol, in which APD 80 and CV were measured as a function of BCL. Fig 11C shows the APDR curve from a simulated monolayer overlaid on data derived from Fig 5 of Bingen et al. [11] Under control conditions the in silico APDR curve corresponds within tolerable limits of the measured data. The simulated curve with blocked I KAChc almost exactly coincides with the experimental curve for the tertiapin-Q-treated cell culture.  [11] for both conditions, in the presence of, and in the absence of the specific Kir3.x-blocker tertiapin-Q. E, AP morphologies at pacing frequencies 1Hz, 2Hz, 4Hz and 10Hz (data obtained from voltage optical mapping of NRAM monolayers). F, AP morphologies at different pacing frequencies, as generated by our monolayer model. The simulated CV restitution curves (Fig 11D) also show good correspondence with the experimentally measured curves derived from Fig 5 of Bingen et al. [11] The changes in AP morphology in response to electrical pacing at gradually decreasing cycle lengths further demonstrate good agreement between experiment ( Fig 11E) and model (Fig 11F).
Spiral wave dynamics in the model and experiments. Next, we studied spiral wave dynamics in our model. The spiral waves were initiated by two independent protocols: (i) the standard S1-S2 cross-field protocol and (ii) by burst-pacing. For measurements of the average period of the spiral, the size of the core, and the trajectory traced by the spiral-wave tip, the spiral was generated in a completely homogeneous, square simulation domain, by the S1-S2 cross-field protocol (Fig 12A). The tip of the spiral wave exhibited complex meandering, resulting in a rounded Z-shaped trajectory with a core-size of 5 mm (Fig 12A, black trace). In monolayers with 17% randomly-distributed myofibroblasts and natural cellular heterogeneity, the tip trajectory retained its shape with local irregularities (Fig 12B). The average period of the spiral wave (as estimated from the dominant frequency of the voltage signal) was found to be 80 ± 4.8 ms (model), which is comparable to 81.3±11.3 ms (experiment) [11]. The average APD 80 during reentry was observed to be 40±4.4 ms (model), which is also comparable to 38.8 ±7.9 ms (experiment) [11]. Wavelength of the spiral-wave, as calculated from the model (CV Ã APD 80 during reentry), was found to be 0.88±0.9 cm, similar to the experimental recording of Fig 12. A, Spiral wave initiated by S1-S2 cross-field protocol in our homogeneous NRAM monolayer model, without myofibroblasts and intercellular variability. The protocol consisted of applying an S1 stimulus of strength 100 pA and duration 2 ms, along the left border of the simulation domain (x 5, 0 y 256), followed by the application of an S2 stimulus of strength 100 pA and duration 2 ms over the region enclosed by 0 x 256 and 162 y 256, when the waveback of the first wave crossed half the spatial extent of the simulation domain. Trajectory of the spiral wave is traced in black. B, Spiral wave in our heterogeneous monolayer model (containing 17% randomly-distributed myofibroblasts and 50-150% intercellular variability) with tip trajectory (black trace). The protocol for spiral-initiation was same as in A. Spiral wave anchored to a clump-type of heterogeneity in the culture dish C and in our monolayer model D, with geometrically identical myofibroblast heterogeneity. E, Phase maps of gradual initiation of reentry upon burst-pacing our in vitro monolayer, with open circle to denote the position of the wavebreak. F, Voltage maps of burst-pacing induced reentry-initiation in our monolayer model with 17% randomly-distributed myofibroblasts and 50-150% intercellular variability.  [11] corresponding to pacing cycle lengths 75-85 ms, which is the range for reentry). In experiments, the tip of the spiral waves were generally non-stationary. They exhibited drift and possibly meander. However, we could not resolve the patterns traced by the tips of the spirals in order to compare with trajectories generated by the model. Therefore, we chose to try and reproduce an experimental situation where the tip of the spiral got anchored to a clump of heterogeneity (Fig 12C). Our monolayer simulations with an identical pattern of myofibroblast heterogeneity also showed anchoring of the spiral wave (Fig 12D). The corresponding wave pattern closely resembled the experimental situation.
To validate the usability of our model in predicting the outcome of in vitro experiments, we burst-paced in silico monolayers at cycle lengths as low as 80 ms and successfully induced reentry. Fig 12E and 12F show in vitro phase maps and representative snapshots of the corresponding in silico voltage maps, respectively, of reentry initiation upon burst-pacing.
To check whether the specific level of introduced intercellular variability, or the presence of electrotonic coupling with myofibroblasts had any influence on the induction of reentry in our monolayers, we performed 2 separate sets of simulations. In the first set, we used 2 different levels of intercellular variability, namely, 75-125%, and 90-110%, alongside the more commonly used 50-150%. At each variability level, we burst-paced 4 configuration-wise different monolayers and observed 100% stable reentry-induction (Fig 13A). In the second set of experiments, we retained 50-150% intercellular variability but removed the electrotonic coupling between cardiomyocytes and myofibroblasts. This ensured that the myofibroblasts acted like inexcitable obstacles. Burst pacing of these monolayers also led to stable reentry induction in 5 out of 5 trials (Fig 13B).
Finally, to validate the usability of the model in predicting results from in vitro experiments, we used our monolayer model to reproduce our experimental previous results from Bingen et al. [11] Previously, we showed that burst pacing of neonatal rat atrial control monolayers led to robust induction of reentrant arrhythmias. Subsequent block of I KACh in these monolayers not only led to the removal of existing spirals, but also prevented spiral-formation by burst pacing [11]. Burst-pacing of our monolayer model also showed arrhythmia induction in 10 out of 10 independent cases (Fig 14A). Once an arrhythmia had been induced, we blocked I KACh completely in our in silico monolayers. Complete block of I KACh resulted in a substantial increase in wavelength through an acute increase in APD. Under the influence of reduced repolarizing currents, the core of the spiral expanded, wavelength increased and the spiral was no longer stable within the confined space of the simulation domain, causing it to disappear ( Fig  14B). Finally, burst pacing of the monolayers with completely blocked I KACh showed that wavebreaks no longer occurred (Fig 14C). The minimum pacing cycle length allowing 1:1 capture of the applied electrical signal shifted to 125 ms, which is in line with the findings of Bingen et al. [11] Therefore, the effect of tertiapin-Q on I KACh-c could also be reproduced with our monolayer model, at different (75-125% and 90-110%) levels of intercellular variability.

Discussion
In the present study, we introduced a mathematical model of neonatal rat atrial tissue, based on available experimental electrophysiological data from NRAMs, and our additional patchclamp data of I KACh in these cells. If data from NRAMs was unavailable, we used parameters suited to NRVMs. We demonstrated the feasibility of our model to successfully reproduce AP characteristics such as dV/dt max , AP amplitude, APD 50 , APD 80 and RMP.
Using our model, we successfully reproduced the APDR and CVR curves generated using optical mapping experiments on monolayers of neonatal rat atrial cardiomyocytes. These curves are essential determinants of stability of arrhythmias in cardiac tissue.
Furthermore, we demonstrated the feasibility of simulations in spatially extended media with our monolayer model by successfully initiating spiral waves by burst pacing. In agreement with the experiments by Bingen et al., [11] our model demonstrated the possibility of removing reentrant circuits upon inhibition of I KACh-c by a specific blocker tertiapin-Q. In addition, our model showed a meandering spiral wave with similar average period (80 ± 4.8 ms, versus 81.3 ± 11.3 ms in experiments), and a rounded Z-shaped core of size 5 mm. Using our monolayer model, we also reproduced spiral wave dynamics from our in vitro data, in the presence of myofibroblast heterogeneity, thus establishing the robustness of the model, for use alongside experiments.
However, in spite of our best efforts to design a state-of-the-art mathematical model, certain inherent flaws remain. A model is only as good as the data that is used to build it. Since the formulations of most of our currents were derived from patch-clamp data from available literature, we faced the problem of dealing with yet incomplete data extracted using different cellisolation, cell-culture conditions and the patch-clamp protocols and procedures for obtaining these crucial data. Regardless of the variability, these factors may influence the experimental outcome to a great extent. Therefore, instead of relying on absolute numbers for the maximal currents from the I-V curves, we used normalized currents to preserve the trend, while compromising in some cases on the absolute values of the currents produced by the cell.
Discrepancies between results obtained from literature and those derived using our model show up in a few important areas. A major limitation of the model is the lack of evidence on Ca 2+ dynamics in atrial cardiomyocytes. For simplicity, we adapted the dynamics for intracellular Ca 2+ from the ventricular model with modifications as described in the S1 Appendix of the supplementary material. However, this lack of key information about Ca 2+ levels in atrial cardiomyocytes can also perhaps account for the slight alteration in model AP morphology. Such difference may impinge on modeling phenomena where accurate representation of the depolarized state is critical. Furthermore, our experimental records of I KACh-c -behavior in single cell versus monolayers raises a fundamental question, that of compatibility between single cell and tissue properties. In this respect, it is important to clarify that being 'constitutively active' is not the same as being 'maximally active'. Even a small, but stimulant-independent activity of a current can be regarded as constitutive activity. In our single cell patch-clamp experiments, we observed a high-conductance stimulant-independent current at hyperpolarizing potentials, which gave way to a low-conductance current at physiological depolarized voltages. However, monolayers constructed using the same cells showed a pronounced response to antagonist tertiapin-Q, which translated into 2-fold prolongation of the action potential, indicating the activity of a high-conductance current in monolayers. In 2006, Cha et al, [62] performed dedicated studies to measure I KACh-c in single and multicellular preparations of canine atrial cardiomyocytes and observed similar anomaly in the behavior of I KACh-c (see Figs 1C, 1D, 2 and 3 of Cha et al. [62]) Our findings are in line with these observations, i.e., modest versus pronounced effect of tertiapin-Q on APD in single and multicellular level, respectively. One approach to quantify I KACh-c involves performing patch-clamp studies on neonatal rat atrial cardiomyocytes in the presence of agonist Carbachol. This technique has been used before for that purpose by other researchers such as Dobrev et al. [9] However, addition of Carbachol to observe substantial difference between the IV curves for basal current and I KACh-c would require (i) an extracellular environment for the cell that contains abnormally high levels of K + , resulting in (ii) a condition that will not reflect the actual situation in monolayers, which we aim to model and study in silico. To this day, this apparent discrepancy between the behavior of I KACh-c at single vs multi-cell level remains unresolved. In fact, this matter may not only involve this one particular current, but is likely to concern a more fundamental aspect of single cell vs tissue vs whole organ behavior regarding ion channel dynamics. Therefore, in this manuscript, we employ a reverse-engineering approach in which we use our single cell patch-clamp data for I KACh-c to build a basal formulation for the I-V characteristic curve in our single cell model. When this model is extended to simulate tissue, we combine the concept of Fig 1A from Dobrev et al. [9] about split-I-V curves at physiologically relevant depolarized voltages and further adjust this formulation on the basis of control and tertiapin-Q-treated APD restitution curves from monolayer studies, to include the strength of the current at varying pacing cycle lengths. Nevertheless, the translation from a single to multi-cell level in terms of ion channel dynamics certainly warrants further investigation. Until new insights become available, we believe, however, that our study introduces a new and valuable model that is robust enough to facilitate simulation studies into key aspects of atrial fibrillation.
A principal limitation of our model is that an experimentally-validated, detailed Q 10 -compensation was not factored in for usability of this model at all temperatures. Because of the limited availability of experimental data on NRAMs, we used a Q 10 compensation similar to Hou et al., [51] with certain other adjustments based on the idea that temperature affects channel kinetics of Na + , K + and Ca 2+ currents as well as their maximal conductances, though to a lesser degree. Because no information was available regarding any of the time constants for the gating variables that controlled the different ionic currents, except fast Na + activation, all the time constants were either taken from the neonatal rat ventricular cardiomyocyte model by Korhonen et al., [41] or adapted to suit the AP characteristic and APDR data fromNRAMs.
Inactivation kinetics of the fast Na + channel is an important determinant of cardiac tissue excitability. In our model, we fitted the steady-state inactivation curve with data derived from Voitychuk et al. [36] The V ½ for this data was~15 mV more positive than in the steady-state inactivation data of Ramos-Mondragón et al. [35] However, we chose to use the data derived from Voitychuk et al., [36] instead of the data derived from Ramos-Mondragón et al., [35] because for a fixed dV/dt max the kinetics from Ramos-Mondragón et al, did not allow wave propagation in a monolayer, as opposed to the kinetics from Voitychuk et al. [36] Taken together, we present in this paper, a quantitatively-robust, computationally-efficient, experimentally-validated mathematical model of the NRAM expressing the I KACh-c . We use this model to build simulation domains that mimic cardiac monolayers in vitro and successfully demonstrate the usability of the monolayer model to study initiation, maintenance and termination of spiral wave arrhythmias by reproducing key results from our previous in vitro experiments.

Methods
Animal protocols were reviewed and approved by the Animal Experiments Committee of the Leiden University Medical Center (LUMC) and conformed to the Guide for the Care and Use of Laboratory Animals as stated by the US National Institute of Health.

Cell isolation and culturing
NRAMs were isolated from hearts of 2-day-old Wistar rat pups as previously described [11,63]. Isolated cells were plated on round glass coverslips (15-mm diameter) coated with fibronectin (Sigma-Aldrich, St. Louis, MO, USA) in 24-well plates (Corning Life Sciences, Amsterdam, the Netherlands). Depending on the assay, cell densities of 0.1-8 x 10 5 cells/well were used and culture as described in our commonly used neonatal rat cardiac myocyte preparation [57,58]. After reaching the giga-ohm seal, the holding potential was set to -50 mV, and the patch membrane was ruptured by gentle suction through the pipette. Whole-cell capacitance (C m ) was calculated from capacitive transient currents evoked during 5 mV steps from a holding potential of -50 mV. The capacitance transients were cancelled with the amplifier. To minimize voltage error and improve the adequacy of the voltage-clamp, pipette series resistance was routinely electrically compensated by >75%. The voltage and current command pulses were applied to the amplifier through its external command input from the A/D converter connected to a personal computer. The instruments were controlled and driven by commercially available Multi-Clamp 700B Commander and Clampex V10.3 softwares (Molecular Devices) for Windows. Throughout experiments, the current and voltage outputs of amplifier were continuously sampled at intervals of 100 μs and recorded onto the personal computer after low-pass filtering at 2-4 kHz with a four-pole Bessel filter.

Solutions
Action potentials and total membrane currents were recorded in the standard extracellular solution containing (in mM): 126 NaCl, 11 glucose, 10 HEPES, 5.4 KCl, 1 MgCl 2 , and 1.8 CaCl 2 (adjusted to pH 7.40 with NaOH). The standard internal pipette solution used for recording of action potentials and total membrane currents contained (in mM): 80 potassium DL-aspartate, 40 KCl, 8 NaCl, 5.5 glucose, 5 HEPES, 5 EGTA, 1 MgCl 2, 4 Mg-ATP, and 0.1 Na 3 -GTP (adjusted to pH 7.20 with KOH). To study the effect of I KACh on AP, 2 μM carbachol (Sigma-Aldrich, St Louis, MO, USA) and 100 nM tertiapin-Q (Alomone Labs, Jerusalem, Israel) were added to standard extracellular solution to activate and block this current in NRAMs, respectively A correction of approximately 11 mV, because of liquid junction potential was applied in the analysis of action potential recordings.
Time-dependent hyperpolarization-activated currents (I KH ) were recorded in a modified extracellular solution. In addition to the standard extracellular solution, tetrodotoxin (30 μM) (Alomonne Labs), nitrendipine (10 μM), 4-Aminopyridine (2mM) and atropine (100 nM) all from Sigma-Aldrich were added in this solution to suppress, I Na , I CaL , I to and muscarinic receptor-activated currents, respectively. The highly selective Kir3.X channel blocker tertiapin-Q (100 nM) was used to isolate I K1 from the total I KH . The standard internal pipette solution used for the experiments contained (in mM): 80 potassium DL-aspartate, 40 KCl, 8 NaCl, 5.5 glucose, 5 HEPES, 5 EGTA, 1 MgCl 2, 4 Mg-ATP, and 0.1 Na 3 -GTP (adjusted to pH 7.20 with KOH)was used to fill the pipettes.

Data analysis and statistics
For off-line analysis, the data stored on the personal computer were analyzed by pClamp V10.3 (Molecular Devices) and GraphPad Prism software version 6 (GraphPad Software, Inc., La Jolla, CA, USA). All data are expressed as means ± standard error of mean (SEM) for a specified number (n) of observations. Significant differences were determined at the P <0.05 level, unless specified. Unpaired Student's t test or the one-way ANOVA test followed by Tukey's test were used for comparing different experimental groups. Statistical significance was expressed as follows: Ã : P<0.05, ÃÃ : P<0.01, ÃÃÃ : P<0.001 NS: not significant.

Construction of virtual monolayers
Virtual monolayers were constructed by allocating cells on a circular selection in a simulation domain containing 256 x 256 grid points. In order to achieve a random distribution of myofibroblasts within the monolayers, N specific sites were selected from within the array, by a random-number generator. For consistency with experiments, N was chosen such that the myofibroblasts constituted 15-20% of the total number of cells comprising the monolayers. The myofibroblasts were coupled to the cardiomyocytes via a gap junctional conductance G gap = 0.5 nS/pF [59]. This coupling is~6 times smaller than myocyte-myocyte coupling. Whole-cell patch-clamp measurements of APs from different samples of NRAMs showed a slight variation in AP characteristics (S1 Fig). When cardiomyocytes as heterogeneous as these are coupled together to form a monolayer, it would therefore be reasonable to expect a certain degree of natural intercellular variability within the monolayer, although most of the heterogeneity in electrophysiological properties would be compensated by electrotonic coupling. Please refer to Fig 7D for proof of variability in the conductance values of currents I K1 and I KACh . It is possible to adopt a very basic approach to accomplish the same in computer models. To incorporate X 1 -X 2 % intercellular variability within the monolayers the following protocol was used: (i) A random-number generator was used to generate 19 numbers from within the range 0.01 x (X 1 -X 2 ), for each cardiomyocyte occupying a site on the cellular grid. (ii) At each cell site, the maximal conductance, associated with an ionic current/flux, was multiplied by one of these 19 random numbers (one random number for each individual ionic current/flux). This ensured a spatially-random distribution between the electrophysiological properties of the individual uncoupled cells. Next, these heterogeneous cells were coupled electrotonically. Electrotonic coupling averages out the electrophysiological properties so that the culture exhibits uniform APD distributions. Simulated monolayers were discarded if the APD 80 maps showed non-uniform distribution of APD or APD dispersion larger than 20 ms. Finally, APD 80 maps from our in vitro cultures were compared with those from the simulated monolayers with varying levels of intercellular heterogeneity. Because the APD 80 maps at 50-150% intercellular variability showed closest resemblance with in vitro data in terms of mean value and spatial distribution, we chose to use this range for our numerical experiments.