A Mathematical Model of the Mouse Atrial Myocyte With Inter-Atrial Electrophysiological Heterogeneity

Biophysically detailed mathematical models of cardiac electrophysiology provide an alternative to experimental approaches for investigating possible ionic mechanisms underlying the genesis of electrical action potentials and their propagation through the heart. The aim of this study was to develop a biophysically detailed mathematical model of the action potentials of mouse atrial myocytes, a popular experimental model for elucidating molecular and cellular mechanisms of arrhythmogenesis. Based on experimental data from isolated mouse atrial cardiomyocytes, a set of mathematical equations for describing the biophysical properties of membrane ion channel currents, intracellular Ca2+ handling, and Ca2+-calmodulin activated protein kinase II and β-adrenergic signaling pathways were developed. Wherever possible, membrane ion channel currents were modeled using Markov chain formalisms, allowing detailed representation of channel kinetics. The model also considered heterogeneous electrophysiological properties between the left and the right atrial cardiomyocytes. The developed model was validated by its ability to reproduce the characteristics of action potentials and Ca2+ transients, matching quantitatively to experimental data. Using the model, the functional roles of four K+ channel currents in atrial action potential were evaluated by channel block simulations, results of which were quantitatively in agreement with existent experimental data. To conclude, this newly developed model of mouse atrial cardiomyocytes provides a powerful tool for investigating possible ion channel mechanisms of atrial electrical activity at the cellular level and can be further used to investigate mechanisms underlying atrial arrhythmogenesis.

Biophysically detailed mathematical models of cardiac electrophysiology provide an alternative to experimental approaches for investigating possible ionic mechanisms underlying the genesis of electrical action potentials and their propagation through the heart. The aim of this study was to develop a biophysically detailed mathematical model of the action potentials of mouse atrial myocytes, a popular experimental model for elucidating molecular and cellular mechanisms of arrhythmogenesis. Based on experimental data from isolated mouse atrial cardiomyocytes, a set of mathematical equations for describing the biophysical properties of membrane ion channel currents, intracellular Ca 2+ handling, and Ca 2+ -calmodulin activated protein kinase II and β-adrenergic signaling pathways were developed. Wherever possible, membrane ion channel currents were modeled using Markov chain formalisms, allowing detailed representation of channel kinetics. The model also considered heterogeneous electrophysiological properties between the left and the right atrial cardiomyocytes. The developed model was validated by its ability to reproduce the characteristics of action potentials and Ca 2+ transients, matching quantitatively to experimental data. Using the model, the functional roles of four K + channel currents in atrial action potential were evaluated by channel block simulations, results of which were quantitatively in agreement with existent experimental data. To conclude, this newly developed model of mouse atrial cardiomyocytes provides a powerful tool for investigating possible ion channel mechanisms of atrial electrical activity at the cellular level and can be further used to investigate mechanisms underlying atrial arrhythmogenesis.

INTRODUCTION
Murine hearts are commonly used as animal models in cardiac research for investigating possible molecular and cellular bases of cardiac arrhythmias (Martin et al., 2012;Huang, 2016). Over years, plethora experimental data on cardiac electrophysiology of the mouse heart have been obtained in a wide range of physiological and pathological conditions. Now it is a challenge to integrate these experimental data obtained at different scales into an integrated mathematical model for systematically elucidating better functional roles of ion channels in atrial electrical excitations and arrhythmogenesis. Up to date, a number of mathematical models have been developed for cellular action potentials (APs) of different types of mouse cardiomyocytes. Examples include the first mathematical model of mouse sinoatrial node by Mangoni et al. (2006), which was further updated by incorporating biophysical properties of membrane ionic currents and intracellular Ca 2+ handling mechanisms by Kharche et al. (2011), the first mouse ventricular cell models developed by Bondarenko et al. (2004), which were updated by incorporating newer experimental data by Li et al. (2010). Recently, mouse ventricular cell models with β-adrenergic signaling regulation and with calmodulin (CaM) mediating Ca 2+ -dependent regulation were developed by Yang and Saucerman (2012) and Morotti et al. (2014), respectively. These models provide useful tools for underpinning insights into the regulation of Ca 2+ handling in physiological and pathological conditions. However, experimental data of mouse atrial electrophysiology is spare due to its limited use for the study of atrial fibrillation, as the mouse heart is too small to initiate and maintain atrial arrhythmias (Huang, 2016) in normal conditions. Recently, however, sustained atrial tachycardia and fibrillation in mouse heart has been reported (Wakimoto et al., 2001), and intensely investigated with transgenic or wild type mice Verheule et al., 2004;Chelu et al., 2009;Choi et al., 2012). These experimental data provide sufficient foundation to develop a mouse atrial cell model. In addition, after decades of experimental studies on the intracellular signaling network (Yang and Saucerman, 2011), sufficient data are also available to construct a biophysically detailed model of mouse atrial myocytes, with consideration of intracellular signaling pathways to simulate the regulation from nervous and endocrine systems (Ramos-Franco et al., 1998;Doi et al., 2005).
Furthermore, electrophysiological heterogeneity within or between the mouse left and right atria has been implicated as a factor initiating and/or maintaining atrial rhythm disturbances (Nygren et al., 2004). Such electrical heterogeneity between left and right mouse atria has also been found in many studies (Lomax et al., 2003a;Trépanier-Boulay et al., 2004;Hu et al., 2005;Nakamura et al., 2010), in which left-to-right differences in several K + currents, as well as in the action potential duration (APD) of isolated atrial myocytes, were reported. A mouse atrial cell model considering the heterogeneity between left and right atria will be helpful on the study of the genesis of AF, such as to investigate how the fast rotors in the left atrium (LA) propagate toward the right atrium (RA) (Sarmast et al., 2003).
The main goal of this study is to develop a mathematical model of mouse atrial myocytes to simulate mouse atrial APs and intracellular Ca 2+ handling mechanisms. In the developed model, Markov chain formalizations were used to describe membrane ion channels for three main potassium currents (the transient outward K + current (I to ), the ultra-rapidly activating delayed rectifier K + current (I Kur ) and the rapid delayed rectifier K + current (I Kr ) and the L-type Ca 2+ channels (LTCCs).
Heterogeneity between left and right atria was considered to derive a left atrial cell model based on the baseline model of right atria. The characteristics of AP and Ca 2+ transient (CaT) generated by our model were compared with mouse atrial experimental findings, which validated the model development. The Ca 2+ -calmodulin activated protein kinase II (CaMKII) and β-adrenergic signaling pathways were also included.

MODEL DEVELOPMENT
Details of model implementation, pacing protocols, and performance of the model can be found in Supplementary Document 1. The source code in C++ of the model is available under the request to the authors (henggui.zhang@manchester.ac.uk).

General Structure
The mouse atrial myocyte model was developed based on the basal model of the mouse ventricular cell model of Morotti et al. (2014) (body temperature), which was extensively modified to match experimental data of mouse atrial cells. The main changes to the original model are described below, and more details are summarized in Supplementary Document 1. The model consists of 159 coupled ordinary differential equations, each of which describes kinetic properties of membrane ion channels, pump and exchanger currents, intracellular ionic homeostasis, and the CaMKII and β-adrenergic signaling pathway. To better describe gating kinetics of ionic channels of I to , I Kur , and I Kr , Markov chain models were used. Equations for CaMKII and β-adrenergic signaling pathways are inherited from the Morotti et al. model, with some parameters of the activation of CaMKII and β-adrenergic activation of protein kinase A (PKA), and their phosphorylation on excitation-contraction coupling (ECC) targets, being modified based on atrial experimental findings. A schematic diagram of the mouse atrial myocyte model is shown in Figure 1.

Modeling of Membrane Currents
LTCC LTCC provides main Ca 2+ influx for myocytes and acts as a trigger for Ca 2+ release from the sarcoplasmic reticulum (SR). There are four types of LTCCs: Ca v1.1 to Ca v1.4 , among which Ca v1.2 plays a major role during cardiac ECC (Hofmann et al., 1999). I CaL has the kinetics of a fast activation and a slower inactivation time course.
The formulation of LTCC in our model was based on a previous Markov chain model developed by Mahajan et al. (2008) (Figure 2A), with CaMKII and PKA phosphorylation modules being updated for mouse ventricles by Morotti et al. (2014). The LTCC model has a minimum scheme of seven states, which incorporates both calcium-dependent inactivation and voltagedependent inactivation, where the transitions between C2 and C1 are strongly voltage-dependent, and the transitions from C1 to O are voltage-independent and determine the steady-state open probability.  2+ and CaM in three compartments. CaMKII is activated by Ca 2+ -CaM binding in the dyadic cleft, sub-sarcolemma area (SL), and bulk cytosol. Phosphatases 1 and phosphatases 2A oppose phosphorylation by either kinase. Ca 2+ and CaM can diffuse across these three regions at different diffusion coefficients. The blue and red circle with "P" in it indicates the phosphorylation of the target by PKA or CaMKII, respectively. I Na , fast Na + current; I NaL , late Na + current; I CaL (or LTCC), L-type Ca 2+ current; I to , transient outward K + current; I Kur , ultra-rapidly activating delayed rectifier K + current; I Kr , rapid delayed rectifier K + current; I ss , non-inactivating steady-state voltage-activated K + current; I K1 , time-independent inwardly rectifying K + current; I KACh , acetylcholine-activated K + current; I KCa , small conductance Ca 2+ -activated K + current; I Nab , background Na + current; I Cab , background Ca 2+ current; I Kb , background K + current; I Clb , background Clcurrent; I ClCa , Ca 2+ -activated Clcurrent; PMCA, plasma membrane Ca 2+ pump; NKA, Na + /K + -ATPase; NCX, Na + /Ca 2+ exchanger.
In order to fit the experimental measures, the voltage dependence of the steady-state transition between C2 and C1 was manually shifted left by 4.4 mV (4.4 mV was added on the voltage in the voltage-dependent transition rates α and β shown in Figure 2A). The current traces shown in Figure 2B were simulated by using the same voltage-clamp protocol as  (Mahajan et al., 2008). (B) Current traces from simulation and experimental records (inset) from Mancarella et al. (2008). (C) Normalized I-V relationship of I CaL . (D) Peak simulated current density compared with experimental data (Rose et al., 1992;Lomax et al., 2003a;Xiao et al., 2006;Mancarella et al., 2008;Xie et al., 2012). All experiments were conducted at room temperature except the data from Xie et al. at body temperature (marked with a red star). described in Mancarella et al. (2008): the I CaL was activated by a series of 250-ms depolarization pulses from a holding potential of −90 mV, to test potentials ranging from −40 to +40 mV (10-mV steps). And the resultant normalized I-V relationship was in good agreements with experimental observations (Mancarella et al., 2008;Hua et al., 2012) (Figure 2C). The total conductance of I CaL was reduced to 1/4 of the ventricular cell model and the simulated peak current density is 4.76 pA/pF, which is in the range of various mouse atrial experimental data shown in Figure 2D.
I to I to is a rapidly activating and inactivating K + current, encoded by K v4.2 and K v4.3 (Nerbonne, 2000). Xu et al. (1999c) reported that the elimination of I to gave rise to substantial prolonged atrial APD, suggesting that I to plays a prominent role in mouse atrial repolarization.
The equations for I to in the mouse atrial myocyte were formulated based on the Markov chain model of Campbell et al. (1993) that comprised seven states as schematically shown in Figure 3A. Transition rates were refitted with a non-manual method (details in Supplementary Document 1). Simulated current traces in Figure 3B were similar to those seen in experiments (the inset in Figure 3B). The current traces were obtained by using a protocol of 500-ms voltage steps in 10-mV increments between −70 and +50 mV from a holding potential of −75 mV. The normalized I-V relationship of our model was then calculated and compared with experimental data ( Figure 3C). Simulation of I to shows a peak current density of 11.6 pA/pF at +30 mV ( Figure 3D).
The curve of steady-state inactivation and recovery from inactivation of I to were shown in Figures 3E,F. To elucidate the steady-state inactivation curve for I to , a two-pulse protocol was applied  (Figures 3E,F). A temperature adjustment factor (Q 10 ) FIGURE 3 | I to in the atrial myocyte model and comparison with experimental data. (A) Schematic diagram of seven-state Markov chain model of I to (Campbell et al., 1993). used to rescale the model to body temperature was determined separately for channel conductance and kinetics. The Q 10 of the channel kinetics were first determined as 1.8 from literature (Courtemanche et al., 1998;Brouillette et al., 2004), and then the Q 10 factor for channel conductance was adjusted and set to 1.1 to match the current density at 32 • C, which is ∼1.34 times of that at 22 • C (Brouillette et al., 2004).
I Kur I Kur encoded by K v1.5 , was identified in atrial myocytes from several species (Lomax et al., 2003a;Wettwer et al., 2004;Pandit et al., 2011). The kinetics of rapid activation and slow inactivation were reflected in the early current nomenclature, such as rapidly activating slowly inactivating current I K,s (Bou-Abboud et al., 2000), or ultra-rapid K + current I Kur (Trépanier-Boulay et al., 2004), which is the currently accepted name for this current.
The HH formulations of I Kur for the mouse ventricular myocyte were replaced by a six-state Markov chain model from Zhou et al. (2012) (Figure 4A), in which the transition rates α and β were decreased to 1/8 of their original values in order to slow down the activation process of the I Kur model to match to experimental data as shown in the Figure 4. The simulated current traces, normalized I-V curve, and peak current density of 5.1 pA/pF at +30 mV were similar to those seen in experiments  Xu et al., 1999c;Brouillette et al., 2004;Trépanier-Boulay et al., 2004;Nakamura et al., 2010;Qin et al., 2012).
(Figures 4B-D) (the current traces was simulated by using a protocol of 500-ms voltage steps in 10-mV increments between −70 and +50 mV from a holding potential of −80 mV). Since all experimental data were obtained under room temperature, the transition rates in the model were fitted to these data then adjust to the physiological temperature. Using the same method as for I to , the Q 10 factor of activation and inactivation transition rates of I Kur was set to 1.9 according to literature (Brouillette et al., 2004). Then the Q 10 factor of channel conductance was chosen as 1.37 to fit the increase in current density at higher temperature (Brouillette et al., 2004;Nakamura et al., 2010), which is about 35 pA/pF at 22 • C and 49 pA/pF at 32 • C found in mouse ventricles.
I Kr I Kr encoded by ERG (KCNH2, K v11.1 ), has been proved to have a key role in the late repolarization of AP in several mammalian species (Sanguinetti and Jurkiewicz, 1990;Wang et al., 1994;Li et al., 1996). The presence and function of I Kr in the mouse atrium has been demonstrated by Nakamura et al. (2010). Another study also shows that the relative mRNA level of ERG channel was higher and the amplitude of I Kr was greater in mouse atria than in ventricles (Mewe et al., 2010).
The HH formulations of I Kr in the mouse ventricular myocyte were replaced with the five-state Markov chain model of Clancy and Rudy (2001) (Figure 5A) with modifications on transition rates to fit experimental measures. The model was stimulated with a voltage-clamp protocol shown in Figure 5B and resultant current traces are in Figure 5C. The time constants for inactivation and activation at several test potentials ( Figure 5D) and the resultant tail current I-V curve ( Figure 5E) were in agreement with experimental data (Nakamura et al., 2010) (the tail current were measured when the clamp voltage returning to −40 mV). The simulated steady-state I-V curve were also shown in Figure 5F. The model was fitted to experiment data at 37 • C, so no temperature scaling factor was needed.
I ss I ss was found in the mouse ventricles and atria by Xu et al. (1999b,c). It is a 4-aminopyridine (4-AP) resistant current, which is distinctively different from I to and I Kur (Brouillette et al., 2004). The activation of I ss is slow and it undergoes little or no inactivation (Xu et al., 1999a).
I ss in our model was reformulated in the form of an HH formula to match its kinetics observed from atrial experimental measures. Since I ss has no inactivating process, one activation gate x ss is sufficient to reflect its relatively slow activation. The simulated current traces and normalized I-V relationship shown in Figures 6A,B were determined by using a protocol of 500-ms voltage steps in 10-mV increments between −80 and +60 mV from a holding potential of −80 mV. Considering the wide range of current density from experiments, conductance of I ss was set to 0.049 nS/pF. The resultant peak current density 3.32 pA/pF at 30 mV and the comparison with experimental measurements are shown in Figure 6C. The simulated peak current density, normalized I-V curve and related current traces are similar to those seen in experiments. Experimental data were mostly acquired at 22 • C, so a Q 10 factor of 2.4 (Brouillette et al., 2004) was applied to the time constant of x ss and channel conductance.

I Kb
In cardiac myocytes, the background potassium current I Kb plays an important role in maintaining the polarized resting membrane potential (RMP) and intracellular K + homeostasis. The recently discovered family of two pore domain K + channels (K2Ps) have properties well suited to a role in mediating background K + conductance (Gurney and Manoury, 2008). The model of I Kb was described by a single Boltzmann equation fitted to experimental data measured by Syeda et al. (2016) at 37 • C ( Figure 6D). They applied 10 mM Ba 2+ in isolated mouse myocytes and attributed the change in background current to twik-related acid-sensitive K + (TASK) channels, which are in the K2P channel family. Since TASK channels display very fast activation and inactivation processes (Gurney and Manoury, 2008), no time constant was incorporated in our model.

I K1
The channel of I K1 encoded by Kir2.x stabilizes the resting potential and is responsible for shaping the initial depolarization and final repolarization of the action potential. The I K1 equation in our model was based on that from DiFrancesco and Noble (1985) with model parameters adapted to mouse experimental data. The 500-ms pulse in the voltageclamp protocol was from the holding potential of −80 mV to voltages from −150 to −40 mV. Currents were normalized to the magnitude at −100 mV. Simulated results were similar to the experimental data (Figures 6E,F).

Other Currents I Na
The I Na accounts for the fast upstroke of APs. Due to the absence of experimental data, in mouse atria, the channel kinetics and conductance were unaltered from the mouse ventricular model of Morotti et al. (2014). The resultant maximal upstroke velocity (dV/dt max ) and AP amplitude were in good agreement with experiments.

I NaL
The kinetics of I NaL were identical to the formulation described in Morotti et al. (2014), due to the lack of experimental data for the mouse atrium. The conductance of I NaL was increased by 68% from the mouse ventricular model to fit the peak current density observed in the experiment (Lemoine et al., 2011).

I KACh
Acetylcholine activates I KACh , an inward rectifying K + current. The molecular composition of cardiac I KACh channels was identified as a heterotetramer consisting of two Kir3.1 (GIRK1) and two Kir3.4 (GIRK4) channel subunits. In the heart, I KACh channels are primarily expressed in the atrioventricular and sinoatrial node and in the atria (Voigt et al., 2014). The I KACh in our model was reparameterized from a HH formulation in a rabbit sinoatrial node model (Zhang et al., 2002). The simulated I-V relationship was compared with mouse atrial experimental data (Supplementary Figure 4).

I KCa
The small conductance Ca 2+ -activated K + current, I KCa , is shown to contribute to repolarization in mouse atrial myocytes (Hancock et al., 2015). Three subtypes of the small conductance Ca 2+ -activated channels (SK1, SK2, SK3) are the carrier of I KCa , which are all sensitive to a selective SK channel blocker apamin (Adelman et al., 2012). We adopted the model of I KCa from another mouse atrial cell model (Asfaw et al., 2020). The channel trafficking effect stimulated by β-adrenergic signaling system was removed.
I NCX NCX (Na + -Ca 2+ exchanger) encoded by the NCX family (NCX1, NCX2 and NCX3), is responsible for maintaining steady intracellular Ca 2+ balance. It extrudes Ca 2+ from the cytosol along with the plasma membrane Ca 2+ pump (PMCA). As a reversible transporter, it also mediates Ca 2+ entry in parallel with LTCCs (Blaustein and Lederer, 1999). Due to the lack of experimental results on NCX electrophysiological data for the mouse atria, the formulations for the NCX were kept unchanged.

I NaK
The I NaK is responsible for maintaining the Na + and K + gradients between the cytosol and extracellular medium, and plays an important role in the restoration of the resting potential, and the generation and propagation of APs (Glitsch, 2001).
The conductance of I NaK was adjusted according to the different regional expression of I NaK in ventricles and atria. Schmidt et al. (1990) found that there were 50% fewer I NaK in the atria than in the ventricles or septa by measuring the number of ouabain binding sites in porcine and canine hearts. In experiments conducted by Wang et al. (1993), the I NaK activity in guinea pig atria was only one third of that in ventricles. In human atrial myocytes, the I NaK activity was nearly 50% lower than that in ventricular myocytes as well (Wang et al., 1996). To adapt this change and also keep the intracellular Na + concentration in physiological range, the maximum rate for the I NaK was decreased by 40% in our model.

Modeling of Intracellular Ca 2+ Handling
The intracellular Ca 2+ handling in this model is based on the framework developed by Shannon et al. (2004)  the SR, the dyadic cleft, the sub-sarcolemma area (SL) and the bulk cytosol. Compared with ventricular cells, atrial cells usually have smaller size, higher percentage of cytosolic compartment, and lower percentage volume of other compartments. Thus, dimensions and parameters for the related structure in our model were modified accordingly (see Supplementary Document 1).

RyR Ca
The formulations for ryanodine receptor (RyR) were based on the model of mouse ventricular myocyte (Morotti et al., 2014), where the half maximal effective concentration (EC 50 ) of the Ca 2+ concentration in SR ([Ca] SR ) on RyR was reduced by 10% to match the fractional Ca 2+ release data from experiments.
The There is study showing that the decay of intracellular CaT in mouse atrial cells is far slower than in ventricular cells (Shanmugam et al., 2011) at 0.5-Hz pacing. Furthermore, Picht et al. (2007) suggested in their experiment in mouse ventricles that the Ca 2+ removal rate via SERCA was greatly elevated at 4 Hz (2.1 times of that at 0.5 Hz) due to increased maximum turnover rate of SERCA (V max ) without change in Ca 2+ affinity of SERCA (K m ). In our model, the V max was reduced by 50% and a coefficient K SERCA describing the effect of CaMKII at high pacing frequencies was adopted which is calculated by: where Ph PLB denotes the level of phosphorylation of phospholamban (PLB) in percentage. And the forward mode Ca 2+ affinity (K mf ) was increased by 13%.

Heterogeneity of Left and Right Atria
Electrophysiological heterogeneity between mouse left and right atria has been investigated in many experimental works (Lomax et al., 2003a,b;Nygren et al., 2004;Verheule et al., 2004;Nakamura et al., 2010). Specifically, the current density differences between left and right atria of I CaL , I to , I Kr , I Kur , and I K1 were measured by Lomax et al. (2003a) and Nakamura There are no significant differences in the recovery from inactivation of I to (Lomax et al., 2003a) and the time constant of activation of I Kr (Nakamura et al., 2010). Considering these experimental findings, a model describing the mouse left atrial AP was derived from the constructed right atrial cell model by increasing the channel conductances of I Kur and I K1 by 110 and 70%, respectively. Validation of the model is shown in the "Result" section. Unless stated otherwise, simulations in this paper were based on the right atrial model.

Validation of the Model
The developed model was validated by comparing the characteristics of simulated action potentials with those of experimental data obtained from mouse atrial cells in several perspectives, including the AP characteristics, CaT characteristics (with their rate dependences). Simulation results on the role of K + currents in the AP were also compared to experimental data for validation or calibration purposes. The rate dependence of intracellular ion concentrations are also shown although cannot be treated as true validation, as the I NaK was calibrated to reproduce physiological [K + ] i and [Na + ] i . Because a large number of experimental data are referred in the "Result" section for validation, in order to avoid confusion, only experimental data explicitly from the mouse atrial cells will be presented in simulation figures for comparison purposes, while the experimental data from other species or other parts of the heart will only be referred in the text.

Simulated AP and Currents
The simulated time course of the mouse atrial AP and major underlying currents are shown in Figure 8A. The AP of left atrial cell model was slightly shorter than that of right atrial cell model due to the significantly greater I Kur and I K1 , while other potassium currents were smaller during the repolarization, producing a partial compensation to the large I Kur and I K1 . Figure 8B shows characteristics of simulated AP under 1-Hz pacing frequency, compared with experimental data. The dV/dt max is about 210 V/s in the atrial model (not shown), which is similar to experimental results (Bao et al., 2016;Syeda et al., 2016). Experimentally measured APDs at 50 and 90% repolarization (APD 50 and APD 90 , respectively) distributed in a wide range, which may be due to different experimental environments, such as temperature, pacing frequency, genetic backgrounds or strains. The computed APD 50 by our model, though was marginally short, matched to the low bound of the experimental data range. The normalized APD based on the ratio of APD 50 over APD 90 were measured and presented (middle panel of Figure 8B), showing that the normalized APD 50 from experiments and our model are in good agreement. Experimental data shown in Figure 8B are listed in Supplementary Table 1.

Intracellular Ca 2+ Handling
Our model well reproduced intracellular CaTs in mouse atrial cells. The CaTs under steady stimuli of 0.5 Hz were depicted in Figure 9A. When paced at 0.5 Hz, the diastolic Ca 2+ concentration in the bulk cytosol ( [Ca] i ) was about 0.11 µM in our model, which was consistent with ventricular experimental findings (Glukhov et al., 2015). As suggested by Walden et al. (2009), there was no apparent difference in the diastolic [Ca] i between murine ventricular and atrial myocytes. The caffeineinduced Ca 2+ transient ( Figure 9B) was similar with experiments by Xie et al. (2012). The effect of caffeine was simulated by increasing the opening probability of RyRs by 7.5 times and complete block of SERCA and I Cab [the same as configurations in the parent mouse ventricular model (Morotti et al., 2014)]. Characteristics of CaTs were calculated from Figures 9A,B and shown in Figure 9C. The peak of amplitude of CaTs were ∼2.6 times higher than the basal level, which is comparable with various experimental data (Li et al., 2005;Shanmugam et al., 2011;Xie et al., 2012). The time constant for the decay phase and time to peak of CaT were also similar with the values reported in Li et al. (2005); Mancarella et al. (2008), and Shanmugam et al. (2011). The fractional SR Ca 2+ release was calculated as the fraction of the maximum [Ca] i before and after application of caffeine.  Figure 9D shows the simulation results of CaTs at various pacing frequencies. Although no direct atrial data available to compare, the amplitude of CaT ( Figure 9E) decreased along with the increasing pacing frequency as seen in ventricular experiments (Picht et al., 2007). When the pacing frequency was between 1 and 3 Hz, there was no obvious change in the maximum [Ca] i , the reduced amplitude of CaT was due to the elevated diastolic [Ca] i . From 3 to 8 Hz, the total SR Ca 2+ release was larger than the Ca 2+ uptake by SERCA, leading to lower [Ca] SR which in turn reduced the amplitude of every Ca 2+ release until new steady states were reached. The new states had a markedly decreased maximum [Ca] i and increased diastolic [Ca] i , therefore, a frequency-dependent attenuation of CaTs.
Ca 2+ fluxes play a key role in the regulation of ECC (Bers, 2000;Eisner et al., 2000). The mouse atrial model can closely resemble the experimental data of Ca 2+ fluxes. Supplementary Figure 2 shows the behavior of RyRs during a cardiac cycle. The peak value of J rel was 95 mmol/(L SR)/s or 2.12 mmol/(L cytosol)/s, which was similar to experimental observation (Lukyanenko et al., 1998). The Ca 2+ influx through LTCC peaked at 0.21 mmol/(L cytosol)/s, which was close to the experimental estimates of 0.30 mM/s (Bers, 2001). The Ca 2+ fluxes of SERCA, NCX, and sarcolemmal Ca 2+ pump are shown in Supplementary Figure 3. In our atrial model, the predominant removal mediator, SERCA, contributed 91.4% of the overall Ca 2+ removal, with NCX significantly dropping to 7.2%. These results were closely similar to the experimental observations in the rat atrium (SERCA and NCX contributed 92.6 and 6.13% of the total Ca 2+ removal) (Walden et al., 2009).

Rate Dependence of Intracellular Ion Concentrations
The rate dependence in concentrations of intracellular Na + , K + is shown Figures 10A,B. For the Na + concentration in bulk cytosol ([Na] i ), although no experimental data for mouse atria, our simulated results well fitted to the values found in murine ventricular myocytes. Without pacing, [Na] i of the model stabilized at about 10.8 mM, which is quite similar to that of 10 mM in the rat ventricle (Despa et al., 2002). It has been found in experiments that [Na] i increases along with the increasing pacing rates. In mouse ventricle, [Na] i was 8.24 ± 4.9, 12.3 ± 4, and 15.1 ± 5.5 mM at 0.1, 0.5, and 3 Hz pacing rates, respectively (Wagner et al., 2006). In rat ventricle, [Na] i was 11.2 ± 2.3 and 15 ± 1 mM at 0.5 and 2 Hz pacing rates, respectively. The [Na] i of our model at various pacing frequencies were well distributed in the physiological range of 10-17 mM. The K + concentration in the bulk cytosol ([K] i ) of our model was about 140 mM at resting, and dropped to around 110 mM when paced at 10 Hz. This change in [K] i significantly influence the equilibrium potential of K + and then the RMP (details in the following section). These results implied that I NaK magnitude was properly implemented, since I NaK is responsible for pumping in K + and Na + out thus keeping the transmembrane gradient of Na + and K + concentration.
The Cl − concentration in the bulk cytosol ([Cl] i ) is also dynamic in our model. It changed from 8 mM to 14 mM when the pacing rate was increased from 1 Hz to 10 Hz ( Figure 10C). The change in [Cl] i generated a higher equilibrium potential of Cl − and played a role in modifying the behavior of Cl − currents (i.e., Ca 2+ -activated Cl − current and background Cl − current in our model). At 1 Hz, there was only an influx of the Cl − forming a transient outward current when the cell was activated. At higher pacing rates, equilibrium potential of Cl − was elevated to about −50 mV, leading to a Cl − efflux in the diastolic phase of AP which slightly depolarized the RMP. Figure 10D shows the rate dependence of the time constant of CaT decay. The decay became more rapid with increasing frequency, which reflected FDAR. At high pacing rates, there was enough CaMKII activated to phosphorylate PLBs (see the section of CaMKII-mediated phosphorylation in Supplementary Document 1 for more details) and enhance SERCA (by a larger coefficient K SERCA described in section "SR Ca 2+ Reuptake (J up )"). This boosted Ca 2+ uptake contributed to the faster drop of CaT.

Rate Dependence of Mouse Atrial APs
We further validated our mouse atrial model against experimental data of AP characteristics at different pacing frequencies (Figure 11). According to experimental data from Bao et al. (2016), the APD adaptation to pacing frequency of mouse atrial cells (shown in Figure 11A) is not distinct, which is quite similar with that in mouse ventricular cells (Wu and Anderson, 2002;Wagner et al., 2009). For restitution, we compared APD 30 after S1 pacing at CL = 130 ms, followed by a single S2 extra-systolic stimulus delivered at various S1S2 intervals (Figure 11B). The APD 30 restitution curve well reproduces the trend of increase when S1S2 intervals were being shortened. The steady-state rate dependence of mouse atrial dV/dt max and RMP are shown in Figures 11C,D. At pacing frequencies lower than 4 Hz (CL > 250 ms), dV/dt max and RMP remained almost stable; at higher pacing frequencies, dV/dt max markedly decreased along with the increase of RMP, which is consistent with experimental data.
There are two likely causes for the elevated RMP at high pacing frequencies: slower late repolarization (defined as the phase of repolarization when membrane potential is below −60 mV) or depletion of [K] i . To investigate which one contributed most to the elevated RMP, we first clamped the [K] i to 140 mM, which is the steady [K] i under 1-Hz pacing, then stimulated the cell with 100-ms CL until the cell reached steady state. In this condition, the final RMP of the cell stayed below −78 mV, which means the rate dependence of the RMP disappeared. Moreover, this relatively low RMP also promoted the recovery of I Na to such an extent that it eliminated the phenomenon of reduced I Na at high pacing frequencies found in experiments by Syeda et al. (2016), and also led to more Na + influx and [Na] i overload. On the other hand, slower late repolarization may also raise the RMP due to not fully repolarized membrane potential when next stimulus comes. The late repolarization was slowed in our model because of reduced I Kr and I ss at high pacing rates. Below −60 mV, I Kr and I ss were reduced by about 50 and 40% at CL = 100 ms compared with those at CL = 500 ms, respectively. According to these changes in current amplitude, we tried to augment the I Kr and I ss to the level of slow pacing (CL = 500 ms) but the resultant RMP (−70 mV) was still elevated. Since it was the clamped [K] i but not faster late repolarization that eliminated the rate dependence of RMP, our simulation suggested the [K] i depletion was the predominant mechanism underlying the increase of RMP when the pacing rate was high.

Rate Dependence of Ion Channel Currents
For better explaining the mechanism underlying the rate dependence of AP, we investigated the role of each main channel current in the model. In Figures 12A,B, ion channel currents under 500, 250, and 100-ms CL were plotted against time or V m . These three CLs were chosen since the rate dependence of currents was unobvious when the CL was larger than 500 ms. There is still no obvious variation of the morphology of the AP with the CL from 500 to 250 ms, whereas distinct changes in I K1 , I NCX , and I NaK could be observed. Due to more APs evoked, much more Na + entered the cell by I Na , resulting in the increased [Na] i and then I NaK . On the other hand, the elevated [Na] i suppressed the I NCX and led to a rise of [Ca] i (as in Figure 9D). I K1 was decreased because its role in the late repolarization was partly taken over by the larger I NaK .
With a further decrease in CL (from 250 to 100 ms), distinct changes on AP morphology could be observed. The RMP was elevated by about 5 mV, resulting from the decrease of [K] i . The overshoot (OS) decreased from 30 to 16 mV, which was mainly due to the decrease of I Na . I Na decrease was also found by O'Hara et al. (2011) in their human ventricular cell model, but the reason of the decrease in our model was the elevated RMP leading to insufficient recovery of I Na and limited Na channel availability, which was different to what was believed to be insufficient time for the I Na to recover. For other currents, because of the smaller OS, I CaL and most currents for repolarization (I to , I Kur , I Kr ) attenuated in varying degrees. It is also worth noting that the decrease in [K + ] i not only elevated the RMP, but also led to a reverse of I K1 due to the altered equilibrium potential of the potassium channel.

Role of K + Currents in the Genesis of Mouse Atrial APs
To further validate the model, the functional roles of main K + currents were investigated by ion channel blocking simulations. Results are compared with experimental data.
Role of I to and I Kur I to and I Kur are both 4-AP-sensitive currents. Studies have demonstrated that relatively low concentrations (≤100 µM) of 4-AP predominantly blocks I Kur while slightly affecting other types of K + currents including I to (Xu et al., 1999b). Since we cannot find experimental data for only I to block, we consider the block of I to and I Kur by 4-AP in a combination manner.
Studies have shown that 4-AP prolongs APs of mouse atrial cell (Lomax et al., 2003a;Nygren et al., 2004;Trépanier-Boulay et al., 2004;Nakamura et al., 2010). Experimental data from 3 different studies are shown in Figures 13Ai,Bi,Ci. Since different amount of 4-AP and pacing frequencies were used in these experiments, to make meaningful comparison with our model, it is necessary to first determine how much potassium currents were affected by different concentration of 4-AP. In experiments, Trépanier-Boulay et al. (2004) found that I to was only blocked 5% by 100-µM 4-AP, whereas Nakamura et al. (2010) found that 100-µM 4-AP may affect I to to a greater extent, so they chose to add 50-µM 4-AP to avoid its effect on I to . We also referred to the dose-dependent data of the effect of 4-AP on I to and I Kur in mouse ventricular cells (Xu et al., 1999b) and assumed that 4-AP should have similar effect in atrial cells. As a result, the percentages of I to and I Kur block were determined to: 16% I to block and 70% I Kur block under 50-µM 4-AP, or 30% I to block and 100% I Kur block under 100-µM 4-AP. The simulated APs and corresponding normalized APDs are shown in Figure 13. The block of I to and I Kur obviously prolongs APD 20 , APD 50 , and APD 90 . The simulation results were in good agreement with experimental data.

Role of I Kr
The role of I Kr in the present model is to modulate the late phase of APs in mouse atrial cells. Normalized APDs and APs with or without application of E-4031 are shown in Figure 14A. I Kr is extremely sensitive to E-4031 and can be fully blocked by only 5 µM of E-4031 (Nakamura et al., 2010). In our simulation, when paced at 2 Hz, the APD 20 , APD 50 , and APD 90 of the cell model were 4.9, 8.1, and 24.4 ms under control condition, respectively. Full block of I Kr didn't change the APD 20 but prolonged APD 50 by 0.3 ms and APD 90 by 2.6 ms. Exposure to E-4031 only affected the repolarizing of APs at levels lower than −25 mV, and it prolonged APD 90 much more than APD 20 and APD 50 . These results were in agreement with experimental findings (Nakamura et al., 2010).

Role of I ss
I ss is not affected by 4-AP lower than 1 mM, but it is sensitive to tetraethylammonium (TEA). Nygren et al. (2004) examined the effect of 5-mM TEA on APDs in isolated mouse atrial preparations. Their results are shown in Figure 14Bi. Since we cannot find dose-dependent measurements in mouse atrial cells, we assume that the sensitivity of I ss on TEA should be similar in the atrium and ventricle. It is reported that 25-mM TEA blocks 58% both I ss and I Kur in mouse ventricular cells (Xu et al., 1999a), but their dose dependences on TEA are quite different. Increasing the concentration of TEA to 125 mM will block I Kur completely but only block 61% I ss in mouse ventricular cells (Xu et al., 1999b), implying that the block of I Kur by TEA has a steeper dosedependent curve and lowering the dose of TEA should significantly decrease the extent of block on I Kur . We did a linear regression on the log scaled dose-dependent curve of the amplitude of I ss and I Kur , and estimated that 5-mM TEA should block about 55% of I ss and nearly no I Kur . Our simulation results were in good agreement with experimental data (Figures 14Bi,Bii).

Novel Achievements
A new biophysically detailed model for the mouse atrial myocyte with considerations of electrical heterogeneity between the left and right atrial cells has been developed, which is the first model representing the spatially heterogeneous mouse atrial APs in the atria. Based on the modification on the mouse ventricular cell model developed by Morotti et al. (2014), the mouse atrial cell model presented here describes cellular Ca 2+ and Na + handling and their regulation by CaMKII and PKA. Three main potassium currents (I to , I Kur , and I Kr ) were developed using Markov chain formalization and fitted to atrial experimental data, which allowed a better representation of atrial specific channel kinetics. Other membrane currents were also updated (I Na , I NaL , I CaL , I ss , I K1 , I KACh , I KCa , I PMCA , I NCX , and I NaK ) for atrial myocytes. All membrane currents were adjusted to physiological temperature either by fitting to data at 37 • C or using temperature adjustment factors. Intracellular Ca 2+ handling (RyRs, SERCA) and the CaMKII phosphorylation module were modified. Considering the heterogeneity between left and right atria, the atrial model was further extended to describe APs of mouse left atrial myocytes. The model was validated by its ability to reproduce the morphology and characteristics of APs and CaTs under various pacing frequencies, including APA, RP, APDs, etc. Characteristics of Ca 2+ transients were comparable with those recorded from mouse atrial cells (see Figure 9). The role of main potassium currents in the genesis of the APs was investigated by simulation and validated with experimental data. The present model provides a basis for further simulating the conduction of excitation waves in the tissue of mouse atria in the future.
In the mouse atrial cell model constructed, we implemented Markov chain models instead of traditional HH models to model some of ion channels. Although the Markov chain model is more complex in the process of parameter fitting, and usually has more differential equations than the HH model and requires more simulation time, it helps to avoid limitations of the HH model in terms of model structure (Rudy and Silva, 2006). On the other hand, the Markov chain model is closely related to the conformation of ion channel proteins, which will be helpful in studying the change of channel function caused by drug action or gene defect in the future.

Comparison With Other Models
The morphology of the simulated AP in our atrial myocyte was relatively short, similar to the report that mouse AP is different from that in most other species (e.g., dog, pig or human) (Kaese and Verheule, 2012). The main factor responsible for short mouse AP is the relative magnitude of the currents present. Mouse myocytes have large outward K + currents, which ensure rapid repolarization. In our model, APD 90 was shorter than that in the ventricular model (Morotti et al., 2014), which was in agreement with experimental observation (Knollmann et al., 2007). In most mammalian species, the atrial AP is characterized by a shorter early repolarization phase, resulting in a triangular wave shape, compared with the ventricular AP with its usually FIGURE 14 | Role of I Kr and I ss in the mouse atrial cell model. (Ai,Aii) Comparison of normalized APDs with experimental data (Ai) and simulated APs (Aii) before and after the application of 5-µM E-4031. The cell was paced at 2 Hz. Experimental data was extracted from Nakamura et al. (2010). (Bi,Bii) Comparison of normalized APDs with experimental data (Bi) and simulated APs (Bii) before and after the application of 5-mM TEA. The CL of the pacing protocol was 130 ms. Experimental data was extracted from Nygren et al. (2004). pronounced plateau phase (Grand et al., 1990). Our atrial model is able to reproduce this triangular atrial AP wave shape. Unlike other mammals, the atrial mouse AP was longer during the early repolarization phase (APD 25 ), compared with that in the ventricle. Our model suggests that it is likely to be the result of reduced current density of I to in atria.
Ca 2+ handling in mouse myocardium has a very rapid mechanism (Bers, 2001). Unlike other mammals with longer APDs, the percentage of Ca 2+ released from intracellular stores (mainly stored in SR) to cytosol during a Ca 2+ transient is much higher in mouse myocytes, and the contribution of Ca 2+ influx via LTCC smaller. Due to the relatively short APD of mouse myocytes, the Ca 2+ removal during diastole is also very fast to ensure a sufficient refilling before the next contraction. Compared with the mouse ventricular cell model (Bondarenko et al., 2004;Morotti et al., 2014), the atrial model presented here has smaller I CaL and I NCX , indicating that the normal ECC coupling requires less Ca 2+ influx and thus less Ca 2+ needs to be extruded from the cell. Therefore, CaTs in mouse atrial myocytes are more dependent on the intracellular Ca 2+ store. Due to rapid SR Ca 2+ release and reuptake, CaTs reach the peak and decay faster than in ventricular myocytes. This is in line with the observation in rat experiments that the time to peak of CaTs is earlier and the decay is faster in atria than in ventricles (Escobar et al., 2004). The difference between atrial and ventricular cells is also reflected in the contribution of three main proteins (SERCA, NCX, and PMCA) to Ca 2+ removal. SERCA contributes more to calcium clearance in the atrial calcium transient, which is also consistent with experimental measurements (Walden et al., 2009).
Recently, another mouse atrial model was developed by Asfaw et al. (2020). Both our model and the Asfaw model provide a set of mathematical equations to simulate the action potential and CaT of mouse atrial myocytes. But there are some distinctive differences between the two as listed below.
(1) Our model considers I KACh , and three main K + currents (I to , I Kur , and I Kr ), which are based on Markov chain formulations, which may be used as a tool for simulating functional impacts of gene mutations and drug interaction in future studies; (2) Our model also considers the electrical heterogeneity of left and right atrial myocytes (Figure 7), and thus extend the model to the left and right atrial cell model; (3) From the perspective of signal pathway, our model includes the effects of CaMKII on a variety of membrane currents, RyR, and SERCA, which plays an important role in the frequency dependent characteristics of the cell; (4) Though both models consider β-adrenergic receptor pathway, while the Asfaw model constructs β1and β2-adrenergic receptors separately; and (5) In the subcellular structure, the Asfaw model considers caveolae and separates the SR into junctional SR and network SR, while our model does not. The simplification of our model on subcellular structure is described in more details in the limitation. From the perspective of model validation, our model was compared with more experimental data. Specifically, we collected the APDs of the left and right atrial myocytes, the characteristics of the CaT, the rate dependence of APD, RMP, dV/dt max , and also the data of the effects of K + current block on the AP. In addition, the normal heart rate of mice is around 600 beats per minute (BPM) (Donner et al., 2011). Our model is stable under 12.5-Hz stimulation (or 750 BPM), while the Asfaw model was only tested with pacing rates up to 6 Hz.
Due to the limited kinetic data on the CaMKII phosphorylation in mouse atrial cells, the CaMKII module is mainly based on the previously published mouse ventricular model (Morotti et al., 2014). Experimental studies have shown that CaMKII can phosphorylate the α1C and the β2a subunits of the LTCC (Grueter et al., 2006). Our model is able to reproduce the quick CaMKII-dependent LTCC facilitation found in many experiments (Valverde et al., 2005;Huke and Bers, 2007). There are also experimental studies showing that CaMKII can phosphorylate the RyR and enhance its sensitivity to Ca 2+ (Wehrens et al., 2004), and the SR Ca 2+ leak and open probability of RyR increase after being phosphorylated by CaMKII (Chelu et al., 2009). Our model is able to mimic these effects by implementing a CaMKII-dependent factor of SR sensitivity and a factor on SR Ca 2+ leak. PLB resides in the SR membrane and negatively regulates the activity of SERCA. PLB phosphorylation will release its inhibition on SERCA, allowing SERCA to more freely uptake Ca 2+ from cytosol. Our simulation shows slow kinetics of PLB phosphorylation (see the section of CaMKII-mediated phosphorylation in Supplementary Document 1 for more details), which is observed experimentally . Compared with the ventricular model, the inhibitory effect of PLB on SERCA in our atrial model is relatively small, because the expression of PLB in murine atrial tissue is only half of that in ventricular tissue (Freestone et al., 2000).

Intracellular Ionic Homeostasis
A valid cell model should be able to not only reproduce the action potential waveform but also ensure the ionic homeostasis at physiological heart rates. Under the prerequisite of keeping 4 types of ions dynamically changing (Na + , K + , Ca 2+ , Cl − ), our model achieves ionic homeostasis under various pacing rates while all ion concentrations were in the physiological range. Dynamically changing ion concentrations enabled us to investigate the relationship between pacing frequency and cell behaviors. For example, we found in simulation that the [K] i depletion at high pacing rates led to a series of influences which was reported by Syeda et al. (2016), including elevated RMP, decreased I Na , and reduced potassium currents. In all the main currents in our model, I NaK and I NCX play a key role of maintaining ionic homeostasis. I NaK is the only current in charge of pumping Na + out and K + in, and I NCX contributes 91% of the efflux of Ca 2+ (the other part is extruded by PMCA). Along with increasing pacing frequency, more Na + influx and K + efflux lead to a higher [Na] i hence larger I NaK until new balance is constructed between the influx and efflux of Na + .
Because of the conservation law, when the steady state is achieved, the net fluxes of each type of ions should be 0 in a whole cycle of an AP. As adopted by other cell models, the stimulus current in our simulation is treated as a K + current to avoid longterm K + depletion. Without doing this, the final balanced [K] i will decrease from 140 to 100 mM, leading to more than 5 mV elevation of the RMP, which is quite prominent. Our atrial cell model was carefully tuned to reproduce reasonable steady-state ion concentrations.
In simulations, a decrease of [K] i was observed during fast pacing. To test whether the rundown of the intracellular K may occur in models as an artifact of the pacing stimulus, we further checked if any change in the stimulus waveform (i.e., different amplitudes of the stimulus current) would affect the intracellular K + concentration. In the model, to evoke an action potential a pacing stimulus current as an analog of K + current with a square wave form, similar to that used in the human ventricular cell model developed by O'Hara et al. (2011) was implemented. In the control case, the amplitude and duration of the stimulus current implemented were 10 pA/pF and 4 ms respectively, which was about 1.5 times of the excitation threshold of the model (at 7 pA/pF with 4-ms duration). In the test conditions, six different amplitudes of stimulus ranging from 80 to 200% of the original stimulus amplitude were implemented at the pacing rate of 10 Hz, and the steady state values of intracellular K concentration were measured. It was shown that up to 200% of the stimulus amplitude used in the study had no obvious effects on the computed [K] i (<5%), suggesting that the rundown of [K] i was not attributed to the artifact of stimulus.

SERCA and FDAR
Frequency dependent acceleration of relaxation (FDAR) in murine hearts has been reported by many researchers (Kassiri et al., 2000;DeSantiago et al., 2002). This phenomenon is related to changes in myofilament properties (Chung et al., 2016) and in CaTs under high pacing rates. The changes of myofilament properties are beyond the aim of this paper, hereby we mainly discuss other possible roles contributing to frequency dependent acceleration of CaT decay (FDAD).
CaMKII plays an important role in the FDAD. The level of activated CaMKII is higher when the pacing is faster, and most targets phosphorylated by CaMKII directly influence the CaT.
Among these targets, SERCA is predominant as it contributes more than 90% in the removal of [Ca] i , i.e., CaT decay, in murine cardiac myocytes (Walden et al., 2009). It is natural to guess that CaMKII amplifies SERCA under high pacing rates through some kind of mechanism. PLB may contribute to FDAD as it is a target phosphorylated by CaMKII and can regulate SERCA. There are reports showing that CaMKII phosphorylates PLB in a frequency-dependent manner (Hagemann et al., 2000;Valverde et al., 2005). But there are also studies suggesting temporal mismatch of the phosphorylation of PLB and the occurrence of FDAD (Valverde et al., 2005;Huke and Bers, 2007), which weakens the view that the frequency-dependent phosphorylation of PLB is the cause of FDAD. On the other hand, Picht et al. (2007) found that FDAD was caused by an increase in the V max of SERCA with presence of CaMKII and unaltered K m (Ca 2+ affinity of SERCA), and was not related to the phosphorylation of PLB. From the perspective of modeling, if we want to guarantee the model has a 200-ms time constant of CaT decay under 0.5-Hz pacing, which is the value found in experiments (Li et al., 2005), the amplitude of SERCA has to be reduced to a degree that is too small at fast pacing rates to pump enough Ca 2+ back in to the SR. In this case, the Ca 2+ concentration in the cytosol and junctional area will be larger, leading to extra RyR Ca 2+ releases during diastolic periods. Considering the exact mechanism of FDAD is still controversy, we chose to amplify the SERCA in an intuitive way that correlates the level of PLB phosphorylation with the V max of SERCA as described in the "Model Development" section. It is proved by the simulation results that our model shows normal FDAD and CaTs.

Lack of Experimental Data
Although the mouse is a widely used animal model in the study of atrial arrhythmias, the quantification of the biophysical properties of mouse electrophysiology is still limited. In some instances, we need to use data obtained either at low temperature (lower than 36-37 • C), from other regions of heart, or from different species, such as mouse ventricles, rats, canines, and rabbits. Most of the ion channels have to be fitted to experimental data acquired at non-physiological temperatures. Those experiments might be conducted at a temperature within the range of 20-37 • C and by different groups. For I ss and I K1 , especially, only data at room temperature were available. The model of I CaL was reformulated to fit with the mouse atrial I-V relationship, whereas the voltage-dependence of inactivation, the time constants of inactivation and recovery from inactivation were not validated due to the lack of available data. Also due to the lack of experimental data, the CaMKII and β-adrenergic modules were adapted from the parent mouse ventricular model. Another issue was the uncertainty about the accuracy of experimental data, which can vary by an order of magnitude, further limited the functionality of our model.

Limitation in the Cellular Structure
In ventricular myocytes, the well-developed transverse tubules (T-tubules) uniformly spread depolarization into the cell, resulting in a nearly synchronous SR Ca 2+ release throughout the entire cell (Cheng et al., 1993). Atrial myocytes in small species, however, either lack a well-developed T-tubule network or have an irregular internal transverse-axial tubular system (Kirk et al., 2003) (Supplementary Figure S1B in Supplementary Document 1). As a result, the initial rise of the Ca 2+ transient starts from the periphery of the cell and then the inner Ca 2+ are released sequentially toward the cell center with decreasing amplitude (Trafford et al., 2013). Since the Ca 2+ handling in our atrial model is based on the framework of ventricular myocytes (Shannon et al., 2004), SR Ca 2+ releases only happen in the junctional area, and Ca 2+ fluxes from the junctional area to SL and SL to cytosol are proportional with the Ca 2+ gradients. This structure cannot mimic the Ca 2+ sparks propagating in the cytosol, therefore underestimates the diffusion speed of Ca 2+ . Increasing the diffusion coefficient cannot be a solution because the Ca 2+ sparks exist in only a short period, which cannot be mimicked by a larger but fixed diffusion coefficient through the whole AP.

Computational Cost
The Markov chain model as a widely used modeling paradigm, its advantages have been discussed extensively (Rudy and Silva, 2006;Fink and Noble, 2009). For example, Markov chain model can better mimic the kinetics of ion channels or proteins, and can be more easily extended to simulate the effects of drugs or gene mutations. However, Markov chain model is generally more complex than HH model, which consists of more differential equations and needs more computation. In the absence of extensive experimental data there are also some limitations in determining the model parameters for the complex system as Markov chain models, and the predictive capabilities of models fit only to experimental I-V curves are limited (Clerx et al., 2019;Whittaker et al., 2020).
Our first intention in constructing the presented model with Markov chain formulations for some ion channels was to describe the electrophysiological properties of mouse atrial myocytes as well as possible, and modulations to them by signaling pathways as detailed as possible, at a cost of non-optimized computational performance of the model. In fact, it costs about 0.8 s for 1 s simulation (pacing the model at 1 Hz) on a 2.8 GHz Core i7 CPU with one core. In multi-dimensional simulation, the computation can be paralleled to significantly reduce the time consumption. The computational time is nearly linear to the number of myocytes, which means that, it may take about 10 s for simulating 1 s electrical activity of a network of 100 myocytes with an 8-core CPU. Since modern CPUs are more and more powerful, we think the performance of our model is acceptable for most applications. If higher performance is needed, one could choose to replace some ion channels with HH models or remove the CaMKII or β-adrenergic pathways from the model.

Conclusion
A new model accounting for the heterogeneous electrical action potentials between the left and the right atrial cells has been developed. The model was derived from, and validated against experimental data obtained from mouse atrial cells. The developed model provides a novel basis for further study of possible mechanisms underlying atrial fibrillations.

DATA AVAILABILITY STATEMENT
All datasets presented in this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
HZ conceived the study. WS, SZ, and HZ designed most of the study, performed the simulations and analyses, and wrote most of the manuscript. WW contributed to the design of simulations, figure design, and manuscript writing. WS aided SZ with reuse of codes from his previous study. HZ, WS, and KW supervised the project. All authors contributed to the article and approved the submitted version.