Human Atrial Cell Models to Analyse Haemodialysis-Related Effects on Cardiac Electrophysiology: Work in Progress

During haemodialysis (HD) sessions, patients undergo alterations in the extracellular environment, mostly concerning plasma electrolyte concentrations, pH, and volume, together with a modification of sympathovagal balance. All these changes affect cardiac electrophysiology, possibly leading to an increased arrhythmic risk. Computational modeling may help to investigate the impact of HD-related changes on atrial electrophysiology. However, many different human atrial action potential (AP) models are currently available, all validated only with the standard electrolyte concentrations used in experiments. Therefore, they may respond in different ways to the same environmental changes. After an overview on how the computational approach has been used in the past to investigate the effect of HD therapy on cardiac electrophysiology, the aim of this work has been to assess the current state of the art in human atrial AP models, with respect to the HD context. All the published human atrial AP models have been considered and tested for electrolytes, volume changes, and different acetylcholine concentrations. Most of them proved to be reliable for single modifications, but all of them showed some drawbacks. Therefore, there is room for a new human atrial AP model, hopefully able to physiologically reproduce all the HD-related effects. At the moment, work is still in progress in this specific field.


Introduction
In the last fifteen years, the increasing interest towards atrial electrophysiology and atrial fibrillation (AF), together with a greater availability of experimental data, led to remarkable developments in human atrial action potential (AP) models [1][2][3][4][5][6]. As a matter of fact, cardiac computational modeling constitutes an efficient tool to investigate the ionic mechanisms involved at cell level and has already been used in a variety of clinical contexts, linking patient manifestations to the underlying electrophysiological mechanisms, thus providing useful insights into different atrial pathologies, including AF, especially whenever experimental measurements were lacking or unavailable [6][7][8][9][10][11][12][13][14][15].
Haemodialysis (HD) therapy represents a unique model to test in vivo, in human, the effects of sudden changes in plasma ionic concentrations and blood volume: in a few hours, patients undergo significant plasma electrolytes variations, together with a significant decrease in extracellular volume. In particular, the HD session causes removal of excess Na + and water, the extent of which depends on the interdialytic weight gain of the patient. Plasmatic K + concentration increases during the interdialytic interval, so that during all HD sessions its level must decrease, while Ca 2+ variations might change depending on the dialysate Ca 2+ concentration and its relationship with pre-HD plasma Ca 2+ levels. [16,17]. These processes often lead to an increased arrhythmic risk for the patient, both during HD and in the hours following the therapy. Indeed, the incidence of AF in end-stage renal disease patients is high: reported rates vary between 7% and 27% [18,19], and HD session may promote AF onset [20,21].
The aim of this paper is first to briefly review the literature concerning applications of the computational approach to the study of the impact of HD therapy on cardiac electrophysiology and then to compare all the currently available human atrial AP models, focusing on their ability to reproduce the electrophysiological changes typically induced by HD sessions, that is, plasma electrolytes and blood volume variations.
In addition, since a modification of the sympathovagal balance in favour of vagal activity may occur during HD sessions in patients showing intradialytic AF episodes [20], the acetylcholine-activated K + current ( KACh ) has been added to all models, and the effect of different acetylcholine concentrations has been considered as well.

Cardiac Cell Modeling and Haemodialysis
Computational models of cardiac AP have already been applied several times to assess the acute effects of HD therapy on cardiomyocyte electrophysiology.
The first attempt in this context was the computational analysis of the heart rate changes during HD [22][23][24]. Since a reliable model of human sinoatrial node (SAN) AP was lacking (as it is still today), these studies were based on a model of rabbit SAN AP, considering the DiFrancesco-Noble model [25,26], as modified by Dokos et al. [27]. Simulation results pointed out that changes of blood K + , Ca 2+ , and pH produce large heart rate variations, showing how electrolyte and pH changes within physiological range may have a remarkable impact on the pace-making rhythm, independently of the autonomic outflow.
The computational approach has been also used to analyse how Ca 2+ and K + changes during HD can alter ventricular repolarization and therefore AP duration [28]. In this work, a model of human ventricular AP was considered [29] and model predictions on AP prolongation were validated against a wide range of experimental data, that is, QT interval prolongation recorded during HD sessions. Simulation results pointed out how computational modeling of ventricular AP may be useful to quantitatively predict the complex dependence of AP duration on simultaneous changes in both Ca 2+ and K + . From this study, a modelbased clinical indication was inferred: Ca 2+ content in the dialysis bath should be designed in order to prevent a critical reduction of serum Ca 2+ , especially in HD sessions with a risk of end-HD hypokalaemia.
The same approach has been applied to atrial electrophysiology: a computational model of human atrial AP has been used to confirm that the intradialytic reduction of plasma K + level is associated with P-wave prolongation [30]. When comparing the simulated atrial APs at the beginning and at the end of multiple HD sessions, imposing in the model the extracellular electrolyte concentrations and heart rate equal to the experimental values measured in vivo, simulation results showed an increase in the time needed for depolarisation and a reduction of the effective refractory period (ERP), both occurring during HD. These two phenomena, in presence of a trigger, that is, repeated premature atrial impulses, frequently induced by a HD session, might form the electrical substrate for intradialytic AF episodes onset. Consistent results were also obtained when performing the same analysis in a multiscale model of the human atrium and considering a simulated ECG [31].
More recently, we applied computational modelling of atrial cellular electrophysiology to the individual case of a patient in which HD regularly induced paroxysmal AF [20]. Simulation results provided evidence of a slower depolarization and a shortened refractory period in pre-AF versus pre-HD conditions, and these effects were enhanced when adding acetylcholine effect in simulation. Starting from these findings, the possible mechanisms leading to intradialytic AF onset were reviewed and reinterpreted. Notably, in a subsequent study, Buiten et al. [21], using the implantable cardioverter defibrillator remote monitoring function, showed that HD is a trigger for AF episodes. In particular, they showed that a lower concentration of K + in the dialysis bath is associated with a higher probability of AF episodes, as predicted by our model-based simulation results.
It is worth noting that, in all these studies, model inputs were set using experimentally measured quantities, that is, plasma electrolyte concentrations and heart rate. However, the actual in vivo extracellular fluid is the interstitial fluid, rather than the blood. Therefore, it could be questioned whether the plasma electrolyte concentrations are a reliable estimate of the interstitial ones, even if this is usually accepted. Indeed, the distribution of free ions between vascular and interstitial compartments has been reported to agree with Donnan theory, which predicts a theoretical ratio between interstitial and plasma concentrations very close to 1 [32].

Computational Models of Human Atrial AP.
Starting from the first two human atrial cell models (Courtemanche et al. [1]; Nygren et al. [2]), both published in 1998, four more have been released in the last few years (Maleckar et al. 2009 [3]; Koivumäki et al. 2011 [4]; Grandi et al. 2011 [5]; Colman et al. 2013 [6]). Hereafter, the six models will be referred to using the initial letter of the first and last authors (i.e., CN, NG, MT, KT, GB, and CZ, resp.). All models consist of a set of ordinary differential equations, each one representing a specific dynamic process occurring in the cell, and the number of equations is related to their complexity: the first models are very simple compared to the most recent ones, where a more detailed description of Ca 2+ handling and cell compartments is included (see Table 1). Moreover, the different parameters and ionic current formulations lead to distinct AP morphologies and properties, for example, AP duration (APD), and CaT duration (CaTD). Since 1998 several papers comparing atrial model performances have been published, mainly concerning CN and NG models, which for many years have been the only ones available [33][34][35][36][37][38][39]. The two most recent reviews [38,39] compared all models except CZ, considering simulations from single cell to whole heart and including both physiological and pathological conditions, thus assessing the current state of the art in atrial computational modeling. Therefore, the comparison of the peculiar properties of these atrial models exceeds the purpose of this work, which rather aims to investigate the acute effects of HD therapy on atrial electrophysiology.
The CN and NG models are almost based on the same human atrial data, and they share most of the transmembrane ionic current formulations: however, CN is developed from the guinea pig ventricular model by Luo and Rudy [40], while NG is developed from the atrial rabbit model by Lindblad et al. [41]. The main differences between the two models are related to Ca 2+ -handling and the CaT is much shorter and with a larger amplitude in NG. As a result, their AP shapes are quite different: a spike-and-dome AP for CN and a more triangular one for NG (see Figure 1, pink and blue traces). The MT and KT models are subsequent extensions of NG: the main changes for MT are new formulations for the transient outward ( to ) and ultra-rapid delayed rectifier ( Kur ) currents, while the KT gives a much more detailed description of Ca 2+ -handling, especially concerning Ca 2+ release. The sarcoplasmic reticulum (SR) is divided into 4 different compartments, including also a spatial dimension: as a result, the CaT is slower compared to the previous model ones, but its duration is increased (see Figure 2, purple trace). The GB model has been developed from the ventricular model published by the same group [42]: most of the ionic current formulations have been preserved and adapted to experimental data acquired in human atrial isolated cardiomyocyte. The AP is quite triangular shaped (see Figure 1, green line), and the Ca 2+ handling is mostly derived from the rabbit ventricular model by Shannon et al. [43], again adapted to human atrial data. It is worth noting that in this model the intracellular K + concentration is kept constant.
The CZ model is the most recent one: it is based on CN, from which he inherited all the ionic current formulations, except for to and Kur , which come from MT. Furthermore, the Ca 2+ handling has been modified using a structure for the SR similar to the one used in KT, together with the corresponding formulation for Ca 2+ release and pumps. Conductance has been slightly tuned, to preserve consistency with the original CN model.
In addition to the models listed above, a different version of the CN model has been considered (from now on referred to as CN * ), slightly modified in order to improve its long term stability [44,45]. This CN * model has been recently used to investigate the specific case study of a HD patient which presented recurrent intradialytic AF, as described above [20].
Moreover, the KT model has been recently modified by the same authors, improving model prediction in chronic AF [14]. The changes involved mostly L-type Ca 2+ current ( CaL ) formulation and this new version of the model (from now on referred to as KT * ) has been considered as well.
Finally, since this study is mainly focused on extracellular electrolyte changes, the known dependence on extracellular K + for both the inward and delayed rectifier K + currents has been added to the atrial models, when not already included [46][47][48][49].
Hereafter, each model will be identified by a specific colour (CN/CN * , pink; NG, blue; MT, cyan; KT/KT * , purple; GB, green; CZ, red) and simulation results for CN and KT will be shown only when a different behaviour with respect to their updated versions (CN * and KT * ) is found.
Simulated APs and CaTs for all the considered models are shown in Figures 1 and 2  Measured CaT amplitudes range from 265 to 345 nM [50,51] and this is best reproduced by CZ and KT. GB produces a slightly smaller CaT while CN, NG, and MT show much higher amplitudes. In addition, the CaT has been reported to decay with a time constant of about 200 ms or even slower [51,53]: such a slow decay is well reproduced by GB and CN only, while in all the other models it is much faster. Model differential equations have been implemented in Matlab (Mathworks Inc.) and a variable order solver has been used to solve them (ode 15 s [54]). Pacing was simulated by a current pulse train (pulses of 3 ms, 1 Hz), maintained for 150 s, in order to allow all the models to reach a proper steady state, that is, intracellular concentrations (Na + , Ca 2+ , and K + ) stable over time. Stimulus current ( stim ) amplitude was set to twice the AP threshold for all models, as previously done in [39] (see Table 1). When using this stimulus, however, the GB model produces an AP which is quite different from the one published in the original GB paper: indeed, some of the biomarkers, for example, AP amplitude and upstroke velocity, are highly stimulus-dependent in this model. Therefore, all simulations with the GB model have been done using the stimulus amplitude needed to preserve the original AP characteristics, which is about 6 times the AP threshold and more close as current density to the ones used for the other models. A summary of the considered models and their main properties is reported in Table 1. At the beginning of the HD session the patient is overhydrated: for this reason, 2-3 litres (or even more) of water is removed from his blood during the treatment. Such a removal is compensated by water refilling from the interstitial fluid and eventually from the intracellular compartment. How fluid accumulation during the interdialytic period and fluid removal during the HD session reflect into variations of intracellular volumes is actually not known in quantitative terms. Therefore, we investigate the effects of a quite large range (±20%) of volume changes. Finally, to explore the effect of vagal stimulation, we added the acetylcholine-activated K + current ( KACh ) to all the models, according to the formulation suggested in [5] and considering the changes induced by 0-15 nM of acetylcholine (ACh).

AP and CaT Biomarkers.
We identified different AP biomarkers in order to quantitatively compare simulation results, choosing in particular the ones already used in previous simulation works either to compare the different atrial AP models [39] or to evaluate the effects induced by electrolyte variations on atrial electrophysiology [30]: action potential duration (APD) was measured as the interval between the AP upstroke and the 90% of repolarization (APD 90 ); resting membrane potential (RMP) was measured at the end of diastole; the AP upstroke duration (AP ud ) was defined as the time needed by membrane voltage to reach 0 mV, starting from the beginning of the pacing pulse [20,30]; AP amplitude (AP amp ) was measured as the difference between the AP peak and RMP; maximum upstroke velocity ( / MAX ) was computed as the maximum slope during the AP upstroke; the effective refractory period (ERP) was measured by simulating a S 1 -S 2 protocol: it has been defined as the longest S 1 -S 2 interval which failed to elicit a S 2 AP of amplitude > 80% of the preceding S 1 AP [55].
A summary of all the considered AP biomarkers is shown in Figure 3, considering the CN * model AP as an example.
In addition, some Ca 2+ transient (CaT) biomarkers have been considered as well, that is, CaT duration (CaTD), measured at 90% of CaT decay (CaTD 90 ), the time needed

Potassium Variations.
In all the considered models the extracellular K + concentration ([K + ] o ) is set to the standard value used in the perfusion bath during V-clamp experiments; that is, 5.4 mM. K + increases during the interdialytic interval and is removed during the HD session: therefore, HD patients often show hyperkalaemia at the beginning of the therapy and hypokalaemia at the end. In order to explore both clinical conditions, we considered the range 3-9 mM. The lower [K + ] o value (3 mM) has been set considering the experimental post-HD measurements available in literature (e.g., 3.9 ± 0.4 [30], 3.6 ± 0.6 mM [56]). The upper [K + ] o value (9 mM) is actually a bit high compared to the pre-HD measurements available (e.g., 4.9 ± 0.5 mM in [30], 5.3 ± 0.9 mM in [56]) but we extended the range on purpose, since there are clinical contexts, such as acute ischemia, in which [K + ] o can locally rise up to 9 mM or more [57].
The main effect of a decrease in [K + ] o should be a hyperpolarization of the cell membrane, due to a different Nernst potential for K + ions. In addition, a [K + ] o decrease leads to a QT interval increase [28], a macroscopic marker of prolonged ventricular APD: therefore, a prolongation of atrial APD is expected as well [58]. On the contrary, ERP is expected to decrease [58,59]: as experimentally observed by Downar et al. [60], APD and ERP may be "uncoupled" when varying [K + ] o . Finally, while slowed cardiac tissue conductivity is a well-known effect of severe hyperkalemia, in the range of [K + ] o concentrations usually measured in HD patients, a positive dependence of conduction velocity on [K + ] o has been observed: this phenomenon is known as supernormal conduction [61][62][63]. Consistently, an increase in PWd during hemodialysis, significantly correlated to K + decrease, has been reported [30]. In a previous simulation study [30] we have shown how both hypo-and hyperkalaemia can cause slowed cardiac tissue conductivity: in hypokalaemia, the RMP is significantly lower (hyperpolarized), and therefore the cell needs more time to reach the membrane potential threshold for AP upstroke; in hyperkalaemia, the RMP is significantly higher (depolarized) and, as a consequence, Na + current availability is decreased and the current is much smaller than usual. In single cell simulations, a slow conduction can be associated to a smaller upstroke velocity and to an increase in the time needed for the voltage to rise toward the AP peak: AP ud and / MAX are then expected to show some kind of U-shape dependence when considering the full [K + ] o range.
A summary of the AP biomarkers for all the different [K + ] o is shown in Figure 4. When some models fail to repolarise with low [K + ] o , the corresponding biomarkers have not been computed.
The models show quite different trends for some of the biomarkers, especially APD 90 and ERP (Figures 4(a)  and 4(b)). In NG, MT, and KT * (Figure 4: blue, cyan, and purple traces), both the RMP and APD 90 behave as expected; however, these models fail to repolarise when [K + ] o is set to low values, exhibiting early after depolarisations (EADs, see, e.g., Figure 5). Indeed, a decrease in [K + ] o leads to a reduction in the conductance of the K + repolarising currents, that is, Kr and K1 , thus prolonging the APD: in these models this effect seems to be overdimensioned, probably due to a low repolarisation reserve, and therefore the membrane potential is not able to go back to its resting value. As an example, in Figure 5, shown are the AP traces corresponding to different [K + ] o levels for the NG model. This is indeed a great limitation when aiming to apply these models to clinical contexts, since normal plasma K + levels are between 3.5 and 5 mM and especially critical for HD patients because they need to remove the K + accumulated during the interdialytic period, primarily in the intracellular pool, and therefore they usually end the HD session in hypokalaemia.
No significant changes have been observed in these models for / MAX and AP ud , while the ERP follows the APD 90 as expected. Finally, the AP amp is inversely related to RMP.
As for the GB model ( Figure 4: green traces), APD 90 and RMP are quite similar to NG, MT, and KT * , but their trends change for low [K + ] o : the model repolarises properly for all [K + ] o , but when considering values lower than 4 mM the RMP is higher (more depolarised) than expected: therefore, APD 90 and AP amp are affected accordingly. AP ud and / MAX show the expected U-shape, related to a reduced conduction for both low and high [K + ] o . As for the ERP, in this model it is always much longer than the corresponding APD 90 , even if it shows a similar dependence on [K + ] o . In addition, probably due to the high stimulus amplitude needed to stimulate the GB model (as explained in the methods section), the ERP could not be computed for most of the [K + ] o : when considering values below 4 mM or above 6 mM, the S 2 AP peak was never lower than 80% of the corresponding S 1 , no matter how short the diastolic interval considered.
As for the CN * and CZ models (Figure 4: pink and red lines), they develop a proper AP for all [K + ] o and they both show a very strong linear dependence of RMP on [K + ] o : this dependence, by itself, should prolong the APD when decreasing [K + ] o , since the membrane potential needs more time to repolarize and then reach its resting value. However, in these models the AP phase 2 shortens as well, so the overall APD 90 is almost constant, or even decreasing with [K + ] o , in contrast with the expected behaviour; this effect is even more pronounced when considering the ERP. As an example, in  In these two models, the AP amp is again inversely related to RMP, and both AP ud and / MAX suggest a reduced conductivity in the [K + ] o range boundaries, especially when considering high [K + ] o .
No significant changes were found in CaT biomarkers and intracellular concentrations, in any of the considered models ( Figure 7).

Calcium Variations.
In all the models the extracellular Ca 2+ concentration ([Ca 2+ ] o ) is set to the standard value used in the perfusion bath during V-clamp experiments, that is, 1.8 mM, which is quite high compared to the normal serum Ca 2+ measured in vivo (1-1.3 mM), as discussed in detail in [64]. During a regular HD session, depending on the dialysis bath concentration, serum Ca 2+ can either rise or decrease. Two previous simulation studies explored the effects of [Ca 2+ ] o on cardiac electrophysiology, considering the range 1-3 mM [64,65]. However, serum Ca 2+ is lower than 1 mM in several patients: reported pre-HD concentrations are, for example, 1.18 ± 0.09 mM in [30] and 1.06 ± 0.16 mM in [28]. Therefore, we decided to extend the explored range to 0.6-3 mM. A summary of the AP biomarkers for all the different [Ca 2+ ] o is shown in Figure 8.
The expected effect of [Ca 2+ ] o increase is a significant decrease of APD [64,66]: the increment in driving force enhances the L-type Ca 2+ current ( CaL ) peak, but at the same time the Ca 2+ -dependent inactivation mechanism is Computational and Mathematical Methods in Medicine  strengthened, thus reducing the overall CaL and therefore shortening the APD. Even if the data showing this inverse relationship between APD and [Ca 2+ ] o has been recorded in ventricular cells, there are a few recordings confirming that the trend is the same for human atrial cells [64]. Indeed, CN, NG MT, KT * , and CZ models are able to reproduce this effect, together with a consistent reduction of ERP (Figures 8(a) and 8(b)). Notably, the original KT model shows an opposite trend for both APD 90 and ERP, fixed in its improved version, where precisely the CaL formulation was changed. On the contrary, GB proves to be not very stable to [Ca 2+ ] o variations: APD 90 and ERP show a biphasic trend, both considerably increasing with [Ca 2+ ] o from 0.6 to 2.5 mM and then decreasing until 3 mM, value in which EADs appear ( Figure 9). RMP, AP amp , and AP ud are almost constant in all models. As for / MAX , only GB, CN, and CZ show a slight linear dependence with [Ca 2+ ] o , related to the increase of CaL peak, which in these models has a greater contribution to the AP phase 0.
As expected [66], diastolic Ca 2+ increases with [Ca 2+ ] o for all the considered models (Figure 10

Sodium and Volume Variations.
In all the considered models the extracellular Na + concentration ([Na + ] o ) is set to the standard value used in the perfusion bath during V-clamp experiments, that is, 130 mM for CN, GB, and CZ, 140 mM for NG, MT, and KT, in agreement with the normal serum levels of 135-145 mM. Na + variation during a regular HD session is usually quite small (e.g., from 139.8 ± 3.4 to 141.6 ± 3.1 in [30], from 129/132 to 133/135 in [20]) and we explored the 120-150 mM range based on the corresponding data available in literature [67,68].
APD 90 and ERP slightly increase with [Na + ] o in all the considered models except GB (Figures 11(a) and 11(b)), while RMP and AP ud are almost constant (Figures 11(d) and 11(f)). In addition, in NG, MT, and KT / MAX increase with [Na + ] o , together with AP peak and therefore AP amp (Figures  11(e) and 11(c)).
No significant differences were found in CaT biomarkers nor intracellular concentrations, apart from an increase in [Na + ] i (not shown). It is worth noting that, in CZ and CN, [Na + ] i regularly shows a stronger sensitivity to changes in extracellular concentrations.
Volume effects have been evaluated by scaling the intracellular volumes of ±20%. The corresponding AP biomarkers variations are all negligible (not shown); for example, KT shows the maximum APD 90 change: +22.5 ms on the whole range. The CaTs become slightly slower when the volume increases, but notable changes have been found only in GB: CaTD 90 increases of +104 ms on the whole volume range, together with an increase of CaT ttp and a reduction of CaT amp (Figures 12(a), 12(b), and 12(c)). In this model [Ca 2+ ] i also increases with volume, while no other significant changes occur in intracellular concentrations (Figures 12(d), 12(e), and 12(f)).
Unfortunately there are no experimental data on the effect of changes in [Na + ] o or volume on cardiac cells to either confirm or deny these findings.

Acetylcholine Effects.
To analyse the effect of a possible increase in vagal activity, we simulated the effects of acetylcholine in the 0-15 nM range, adding to all models the same KACh formulation used in [5]. The expected effect of an additional outward K + current is a more hyperpolarized RMP, together with a shortening of APD and ERP. This has been confirmed by experimental data [69][70][71] as well as by previous modeling studies [5,72,73]. When considering concentrations higher than 3 nM, all the considered models show a significant decrease of both APD 90 and ERP (Figures 13(a) and 13(b)). In addition, the resting membrane potential is indeed hyperpolarised, especially in NG and MT (Figure 13(d)). The AP ud is inversely related to RMP changes ( Figure 13(f)), whereas AP amp and / MAX keep almost constant (Figures 13(c) and 13(e)); therefore, the overall conductivity is slowed down by ACh. In fact, starting from a more hyperpolarized potential with no significant changes in / MAX , the cell needs more time to reach the threshold for Na activation, to produce the upstroke. On the contrary, in GB both AP amp and / MAX increase with ACh, mostly due to a larger Na + current for lower RMP, thus compensating this effect and limiting the theoretical AP ud increase and the corresponding reduced conductivity.
[Na] i (mM) [K] i (mM) Negligible effects were found in CaT biomarkers and intracellular concentrations (not shown).

Discussion and Conclusions
We have briefly pointed out that computational models of cardiac action potential (AP) have been successfully applied to investigate HD-related effects on the electrophysiology of different cardiac tissues (sinoatrial node, ventricle, and atrium) often leading to relevant interpretations of macroscopic observations made in clinical ECG and/or useful suggestions about HD treatment personalisation.
However, all these studies have been performed by using cardiac cell models that had been developed on the basis of in vitro experimental data, almost always acquired using standard Tyrode's solutions as extracellular fluid. It is obviously correct to simulate the electrical activity of cardiac cells by imposing the same conditions used in experimental protocols as far as the aim is a comparison with in vitro experimental data. On the contrary, it can be incorrect to use the same constant concentrations when the ultimate aim of simulations is the analysis of in vivo, dynamical conditions such as those during a HD session. Sometimes, this possible cause of discrepancy has been mitigated by few changes to the original models, for example, introduction of the effect of extracellular   pH on the Na + /K + pump activity in the DiFrancesco-Noble model of SAN cell [22] and strengthening of the CaL Ca 2+dependent inactivation in the Ten Tusscher-Panfilov model of human ventricular cell [28].
However, a systematic analysis of the applicability of cardiac cell models to reproduce the specific conditions occurring during HD or, in general, when the extracellular fluid composition changes is still lacking.
In the present paper, we addressed this kind of problem by focusing on human atrial cell models and on the following "cell environment" changes: extracellular electrolyte concentrations (K + , Ca 2+ , and Na + ), cell volume, and acetylcholine.
We pointed out that several human atrial models are available, with significantly different behaviour upon such environment changes.
Unfortunately, experimental data on human atrial cells induced by extracellular concentrations changes are really rare in literature. This makes a stringent quantitative comparison between simulations and experimental measurements not possible for most of the considered electrophysiological properties. On the other hand, some qualitative behaviour is expected based on the overall evaluation of (i) knowledge of physiological mechanisms (e.g., the link between membrane resting potential and Nerst K + potential); (ii) in vitro data measured in different cell types and species (e.g., [50,64]); (iii) in vivo data on macroscopic ECG markers known to be related to atrial cellular electrophysiology (e.g., PWd). We found a major problem in the NG, MT, and KT models: they all fail to repolarize and to produce physiological APs when [K + ] o is lower than 4 mM. This makes these models not appropriate to simulate the cardiac impact of HD. Indeed, the change in plasma [K + ] o is one of the more important and quantitatively large effects of HD, since K + removal is one of the treatment aims and the end-HD [K + ] o is almost always much lower than 4 mM [30,56]. Indeed, even in control condition ([K + ] o = 5.4 mM), the repolarising K + currents (especially Kr and Ks ) of these models are quite tiny when compared to the ones of CN * or CZ, who repolarise properly up to [K + ] o = 3 mM (current peaks are about 10 times smaller). In the NG paper the authors explicitly say that the Kr conductance has been reduced to fit AP data and that this current has been assigned a very low density [2]. Therefore, an increase of its conductance may improve the performance of these models for low The GB model exhibits several shortcomings as well. First of all, although it produces a proper AP at all the tested [K + ] o , it behaves nonphysiologically when [K + ] o is lower than 4 mM: the RMP depolarizes instead of hyperpolarize and, as a consequence, the APD 90 also goes in the opposite way (shortening) and / MAX dramatically decreases. Moreover, the excessive sensitivity to the amplitude of the stimulus current makes the computation of the ERP very unstable, leading to too long ERP values or no ERP at all. Since also in GB the 14 Computational and Mathematical Methods in Medicine repolarising K + currents (both Kr  ] o is lower than 1 mM. Therefore, the GB model turns out to be completely unsuitable to simulate the HD conditions. As for the wrong dependency of APD on [Ca 2+ ] o a possible solution should address a modification of the L-type Ca 2+ current, increasing the Ca 2+ -dependent inactivation with respect to the voltage-dependent one, in order to reduce the overall current for higher Ca 2+ levels, despite the increase in driving force. Indeed, Ca 2+ -dependent inactivation seems to be underestimated in many AP models [64], and previous modeling works managed to reproduce the inverse APD-[Ca 2+ ] o just by strengthening this mechanism [65,74]. The CN * model responds properly to [K + ] o changes, at least from a qualitative point of view. It also reproduces well the "uncoupling" between APD and ERP variations when [K + ] o is increased (APD slightly decreases whereas ERP increases): this was experimentally reported by Downar et al. [60] when perfusing cardiac cells with hyperkalaemic RMP: resting membrane potential; APD: AP duration; ERP: effective refractory period; APud: upstroke delay, inversely correlated with conduction velocity; black arrows: expected increase/decrease during HD, with the corresponding references, +/++/−/−−: moderate/large biomarker increase/decrease; * : biomarker changes different from the expected ones.
"ischemic blood" and interpreted as a secondary effect to changes in resting potential, which is known to affect, in turn, the Na + channels. In addition, simulation results for the CN * model predict a decrease in intracellular Na + when increasing [Ca 2+ ] o : we are not aware of any available experimental data to confirm/deny this observation, which could have relevant implications.
The CZ model exhibits a good stability, with neither repolarization failure nor EADs occurrence. However, it also has a few discrepancies with respect to the expected behaviour: APD decreases with [K + ] o , while the opposite should happen [28] and CaT amp is insensitive to [Ca 2+ ] o while in all the other models it increases with it, in agreement with experimental data reported in [66].
As for the quantification of cardiac side-effects of HD therapy, overall simulation results confirm that changes in [K + ] o and [Ca 2+ ] o are the ones mostly affecting cellular electrophysiology [28,30], whereas [Na + ] o and volume seem to have a minor impact. A qualitative summary of the expected variations in [K + ] o and [Ca 2+ ] o during HD and of the corresponding biomarker changes is shown in Tables 2  and 3, respectively, comparing experimental/computational data from the literature with the simulations results of this study.
Simulation results of acetylcholine effect show a reduction of APD and ERP in all the models, together with a more hyperpolarised RMP, in agreement with experimental data and previous modeling studies [5,[69][70][71][72][73]. In addition, all the models except GB show a reduction in AP ud , suggesting a slower conductivity, also consistent with the increased vulnerability to arrhythmias, such as AF, due to an increased vagal activity [20,70]. However, there is no experimental evidence to confirm or deny these results, and a more detailed description of autonomic regulation should be considered for future improvements in computational modelling of acetylcholine effects.
Other HD-related effects (e.g., acidosis correction) have not been addressed in our analysis and are left to further investigations.  Finally, it is worth remembering that HD patients are first of all uremic patients: this pathological condition (e.g., "uremic intoxication") can also affect some aspects of cardiac cellular electrophysiology and should be incorporated into the models. As a relevant example, downregulation of the Na + /K + pump and high levels of circulating Na + pump inhibitor have been reported in uremic patients compared to individuals with normal renal function, by several investigators [75][76][77][78][79][80].
In conclusion, computational modeling of human atrial cells constitutes a very useful tool to investigate the electrophysiological changes occurring in patients undergoing HD therapy. Nevertheless, it is always important to select carefully the specific model to use, depending on the particular aspect of interest. Currently, CN * seems to be the more suitable human atrial model to analyse HD-related effects on atrial electrophysiology, though it is the oldest one and, therefore, it has a less detailed description of several cellular mechanisms.
Therefore, an additional model could be developed, trying to integrate and reconcile the knowledge of cellular and subcellular processes and their reactions to changes in the extracellular environment, taking into account the possible suggestions given above.
In this respect, work is still in progress in this specific field.