A Computational Model of Electrophysiological Properties of Cardiomyocytes

Introduction. The method of electrical analogies for the analysis of bioelectric dynamic processes in cardiomyocytes is used in the study. This method allows for replacing investigation of phenomena in non-electrical systems by research of analogous phenomena in electrical circuits. The investigation of time processes in cardiac cells is based on the solution of the system of ordinary differential equations for an electrical circuit. Electrophysiological properties of cardiomyocytes such as refractory period, maximum capture rate and electrical restitution are studied. Mathematical modeling. Computational simulation of the action potential and currents for 𝐾 + , 𝑁𝑎 + , 𝐶𝑎 2+ ions in cardiomyocytes is performed by using the parallel conductance model. This model is based on the assumption of the presence of independent ion channels for 𝐾 + , 𝑁𝑎 + , 𝐶𝑎 2+ ions, as well as leakage through the membrane of cardiac cell. Each branch of the electrical circuit reflects the contribution of one type of ions to total membrane current. Results.The obtained electrical restitution curves for ventricular and atrial cardiomyocytes are presented in the paper. The proposed model makes it possible to identify the areas with the maximum slope on the restitution curves, which are crucial in the development of cardiac arrhythmias. Dependences of calcium current on stimulation frequency for atrial and ventricular cardiomyocytes are obtained. Analysis of the kinetics of calcium ions under various protocols of external influences can be useful for predicting the contractile force of cardiomyocytes. Conclusion.The results of calculations can be used to interpret the experimental results obtained in investigations of cardiomyocytes using the “laboratory on a chip” technology, as well as in the design of new experiments with cardiomyocytes for drug screening, cell therapy and personalized studies of heart diseases.

Conclusion.The results of calculations can be used to interpret the experimental results obtained in

Introduction
The method of dynamic analogies is widely used for a long time as a basis for interdisciplinary research of technical systems and physical models in medicine, biology, ecology [1][2][3][4]. The concept of biomimetic design [5] is the development of dynamic analogies methods.
In electrical circuits, electrical energy is transmitted through the branches containing resistors, capacitors, inductive coils and other components, and redistributed between branches by means of nodes. Electrical processes, including nonlinear ones, are investigated using known concepts: electric current, voltage, electromotive force. The mathematical description of electrical processes often coincides with the description of processes in objects of a different physical nature, which allows us to replace the study of phenomena in non-electrical systems by studies of analogous phenomena in electrical circuits. Analyzing the components and topological equations of various types of systems, we can detect their dynamic analogies. This makes it possible to use Kirchhoff's laws of electrical engineering and also component equations, in particular, for analyzing bioelectric processes in living tissues and cells. In this case, the cell membrane is represented by a circuit model that includes a capacitive element, and in which the ion channel conductivities are represented in the form of resistive linear and nonlinear components, but nonequilibrium electrochemical processes are described by voltage sources. This approach was applied by Hodgkin-Huxly during the study of bioelectric processes in the nervous tissue [6].
It should be noted that Kirchhoff's laws are basic for analysis of electrical model that reflects bioelectric processes in cells. Application of the Kirchhoff's laws and component equations to electrical circuit leads to system of differential equations. Thus, the study of action potentials and time processes in cardiomyocytes can be based on the system of ordinary differenti-al equations for a nonlinear model of an electrical circuit that describes the dynamics of cardiac cells functioning.
1 Literature review and problem statement Nowadays, lab-on-chip technology is an emerging in vitro tool to obtain and study cardiomyocytes typically with human-induced pluripotent stem cells (hiPSCs). These cells can be differentiated into a variety of cardiomyocytes (hiPSC-CMs) and then used for the development of heart disease models, drug screening and tissue regeneration for cell-replacement therapies [7][8][9][10]. The key conclusion of these studies is the similarity of electrophysiological properties in hiPSC-CMs and human cardiomyocytes.
However, experimental studies using hiPSC-CMs are fairly complex and have advantages, as well as limitations, these latters being mainly related to the maintenance of the environmental and stimulation conditions through time, the small numbers of myocytes produced and their immature nature. As a result information about electrophysiological properties of hiPSC-CMs is limited.
In this scenario, further improvements of methods and tools to study hiPSC-CMs' electrical activity can be provided by computational modeling of electrophysiological properties of cardiomyocytes.
In order to investigate the functional properties of hiPSC-CMs, different electrophysiological technologies are used. Studies [11,12] are based on patch-clamp technique, which is a classical approach to record intracellular electrical activity by inserting a sharp electrode into a cardiomyocyte. Alternatively, multielectrode potential recordings [10,12] allows for studying the ion channels functions. Contactless imaging methods with use of voltage-and calcium-sensitive fluorescent dyes [13,14] provide multicellular recordings at high spatial resolution for studying morphology of transmembrane potential and intracellular calcium transients during cardiomyocyte's differentiation and drug discovery. The non-invasive, high-resolution method [15], based on genetically encoded calcium and voltage fluorescent reporters, allows for the longterm study of healthy or diseased hiPSC-CMs and, сonsequently, mechanisms of arrhythmia.
Noteworthy, the study of hiPSC-CM electrophysiological properties allows for the assessment of the functional maturity of the cells. In agreement with [11,13], three action potential types (nodal-, atrial-, or ventricular-like) are generally identified.
A variety of bioengineering strategies, employed with different cell types, can prove the influence of different factors on hESC-CMs functionality. In [16,17] the authors have recorded membrane potential of cardiomyocytes and determinated the action potenti-al's duration after electrical stimulation and under spontaneous beating. In [18] the design and fabrication of a micro-scale cell stimulator, capable of simultaneously providing mechanical, electrical and biochemical stimulation, have been described.
To recapitulate physiological environment of cells in the native myocardium we have recently developed specific heart-on-chip technologies [19]. The proposed device allows for performing electrical and mechanical stimulation and for evaluating the electrophysiological properties of cardiomyocytes. Major attention is paid to the study of the functionality of micro-cardiac tissue: the electrical functionality of control and stimulated constructs was assessed in culture by evaluating the excitation threshold and the maximum capture rate.
However, most of the quoted works describe the characteristics of cardiomyocytes at the sarcomere level, myofilament organization, ion channel expression and intercellular connections during the differentiation process.
The present study is devoted to computational modeling of cardiomyocytes' electrical activity that supplements experimental modeling and can predict the kinetics of contractile force of heart cells at the functional level.

The aim and objectives of the study
The aim of the paper is to focus on the cardiomyocyte electrophysiological properties at the functional level, including the generation of action potentials, activativation/inactivativation processes in calcium ions channels, the frequency-dependent changes in action potential duration and the intracellular calcium release or uptake, that enables to explain the changes in excitation-contraction coupling of cardiomyocytes.
The experimental modeling of electrical and mechanical processes in cardiac cells, realized on the lab-on-chip platform [19], was a starting point for the development of the purely computational model.
Considering the peculiarities of the research on the heart-on-chip platform, our work was focused on: the improvement of the parameters of the mathematical model of cardiomyocytes for reflection order to recapitulate the functional maturity of the heart cells; the study of the cardiomyocytes' refractoriness phenomenon; the investigation of the processes involved in the cardiomyocytes' electrical restitution.

Computational modeling
In this paper, the parallel conductance model was used to modeling of cardiomyocytes' electrophysiological properties at the functional level. This model, based on the approach of Hodgkin-Huxly for nerve tissue, was improved by many scientists specifically to address cardiomyocytes [20,21].
The proposed model [22] accounts for main ionic currents and change of potential for cell's membrane. Independent conductance channels are used for + , + , 2+ and leakage. Transmembrane potential ( ) can be represented as the sum of alternating component of membrane potential ( ) and resting potential 0 : As it follows from the parallel conductance model, the alternating component of the membrane potential is described as: ≥ is the depolarizing pulse of current (with amplitude 0 and duration ), which is repeated with a stimulation frequency ( ), are currents for potassium, sodium and calcium, respectively, is the resting potential, 0 , 0 and 0 are the conductances of potassium, sodium and calcium ions at rest; , and are Nernst potentials of potassium, sodium and calcium respectively; is the leakage conductivity through the membrane; is the electromotive force of source, which simulates Nernst potential for chlorine ions, leakage and other factors that affect the membrane potential at rest.
Membrane conductances for + , + , 2+ channels are described by the following equations: where max , max and max are conductances for potassium, sodium and calcium ions, respectively, in the case that all the channels for this type of ions are in the open state; is activation function of + channels; is activation function and ℎ is inactivation function for + channels; is activation function and is inactivation function for 2+ channels. Consequently, conductances and currents for + , + , 2+ ions are determined by five gating variables , , ℎ, , and , which are solutions of the differential equations: where ∞ , ∞ and ∞ are the steady-state value of activation function for potassium, sodium and calcium channels respectively; ℎ ∞ and ∞ are the steady-state values of inactivation function for sodium and calcium channels; , and are the relaxation periods of activation for potassium, sodium and calcium channels; ℎ and are the relaxation periods of inactivation for sodium and calcium channels.
Equations (1) and (3) define the Cauchy problem for the system of ordinary differential equations with the initial conditions: where 0 , 0 , 0 , ℎ 0 and 0 are the steady-state values of activation and inactivation function if alternating component of membrane potential is equal to zero.
The attained system is a set of stiff differential equations. Consequently, to solve Cauchy problem, the implicit methods of integration is used [23].
Using the predictor-corrector method, step ∆ in the initial segment of integration should not be too large (∆ < , where is duration of the depolarizing pulse).
A detailed description of the functions and numerical values for parameters of the proposed model is given in [22].

Numerical experiments
Numerical experiments to model AP in heart cells were performed in Matlab environment. The generation of AP and ion currents for ventricular and atrial cardiomyocytes was simulated. The study accounts for the cardiomyocyte electrical properties such as the refractory period, the maximum capture rate and the restitution during stimulation.

4.1
Action potentials and main currents for cardiomyocytes Simulated action potentials for ventricular and atrial cardiomyocytes (Fig. 1) were obtained in [22]. In accordance with the theoretical information [20] each action potential has three characteristic phases: depolarization, plateau and repolarization. Depolarization phase is determined by a sharp increase of AP amplitude initiated by the growth of membrane permeability for sodium ions. Plateau phase describes the slow decline of action potential, because the inward (slow calcium) and outward (potassium) currents are nearly balanced (Fig. 2a). Repolarization phase is characterized by the faster decline of AP, which could be explained by inactivating of calcium channels and activating of the potassium channels [20].
According to the proposed model [22], currents and conductances for potassium, sodium, calcium ions were calculated and presented for ventricular cardiomyocytes in Fig. 2b-d (currents ions, conductance of 2+ channels and conductance of + , + channels, respectively). In numerical experiments the stimulation frequency (stimulation cycle length ( )) varied from 1 Hz (1000 ms) to 6 Hz (167 ms). Durations of action potentials (APD) were measured as the interval from beginning of action potential to the time point of 90% of AP duration (APD90) (Fig. 1).

4.2
Refractory period and maximum capture rate of cardiomyocytes It is known [20] that refractoriness of cardiomyocytes is determined by the refractory period (RP), which is a period of time when another AP cannot be generated as the cell needs to recover from the previous AP. This peculiarity of AP takes place due to the inactivation of + and 2+ channels. Three action potentials are chosen as representative APs at various stimulation frequencies (2.5 Hz, 3 Hz, 3.25 Hz, 3.27 Hz). Fig. 3a Different authors have studied the electrical restitution that is an intrinsic heart property of change of action potential duration (APD) according to heart rate (HR) [24][25][26][27][28]. The relationships between APD of cardiomyocytes and (or ) or, more correctly, preceding diastolic intervals ( ) are investigated. Usually, the preceding is determined as the time interval between full cardiac repolarization and the growth of the next AP (during cycle length).
The APD shortens with decreasing length and thus with decreasing . It is thought that restitution takes place because calcium current does not fully recover at short , which leads to short APD at short . Electrical restitution data were simulated using the dynamic restitution protocol (DYRT) [24], in which stimulation impulses were generated with various stimulation frequencies or cycle lengths . According to this protocol cardiomyocytes were stimulated in the physiologically determined frequency range with increasing incrementally until the myocyte fails to capture AP, that means refractory period has been reached and maximum capture rate of cardiomyocytes has been obtained.
The responses (AP and currents) from ventricular cardiomyocytes stimulated at different were investigated. Electrical stimulation started at 1 Hz and increased up to 6 Hz with assessment of the maximum capture rate. The ventricular myocyte failed to capture AP at stimulation frequency = 3.27 Hz ( =306 ms).
Similar action potentials and currents in atrial cardyomyocytes were obtained varying . However, the atrial myocyte failed to capture AP at higher stimulation frequency compared to the stimulation frequency of ventricular myocyte (at frequency =5.125 Hz, when maximum capture rate is achieved).
APDs for atrial and ventricular cardyomyocytes shortened at increasing stimulation frequencies and decreasing cycle lengths (Fig. 4). Mean APD90 at each stimulation frequency (cycle length) was obtained by averaging of 20 consecutive steady-state APs.
The restitution curves of mean APD90 for atrial and ventricular cardyomyocytes are shown as functions of stimulation frequency (Fig. 4a) and as functions of cycle length, the inverse of stimulation frequency (Fig. 4b).
Сalcium current mean value for atrial and ventricular cardyomyocytes decreases with increasing of stimulation frequency (Fig. 5a) and with decreasing of cycle length (Fig. 5b). Mean value of at each stimulation frequency (cycle length) was obtained by averaging of during 20 consecutive action potentials.

Discussion
Many studies have been performed in order to investigate the changes of APDs, which likely plays an important role in arrhythmogenesis. According to the overview of published data [24][25][26][27][28], the normal electrical restitution curve has three main phases: a steep recovery at the shortest either or , a transient decline at middle ( ), and a final rise to a plateau at long ( ). Furthermore, the maximum slope of the restitution curve is crucial in determining the arrhythmogenic properties of cardiomyocytes. Under maximum slope of the curve more than 1 the fluctuation of cardiomyocyte's electrical activity occur thus creating the heart unstability. In addition, rate dependent alterations of APD are markers of atrial and ventricular arrhythmias.
Dependences APD against (Fig. 4a) and APD against (Fig. 4b) generate the APD electrical restitution curves, which also have several phases with the various steepness. The maximum slope of the curves arises due to refractoriness of cardiomyocytes and corresponds to the range of maximum capture rate. Therefore the analysis of the electrical restitution curves allows for explaining the difference between behavior of atrial and ventricular cardiomyocytes.
The results of the numerical experiments were analyzed and compared with the experimental results, based on the lab-on-chip technology [19]. In the computational analysis the same research protocol was used as for the experimental study. According to the   dynamic restitution protocol, cardiomyocytes were stimulated at the physiologically determined frequency range with increasing incrementally until the myocyte fails to capture AP. This means that the refractory period has been reached and the maximum capture rate of cardiomyocytes has been obtained. During experimental results, the values of maximum capture rate were obtained in the frequency range from 3 Hz to 5.5 Hz. This variability can be explained by the presence of different types of cardiomyocytes (atrial-or ventricular-like) after the differentiation process from human-induced pluripotent stem cells.
The calculated restitution curves allows for identifying the maximum slopes, which determine the arrhythmogenic properties of heart cells (in which cardiomyocytes can change the beating rate). Moreover, there is supportive evidence that drugs, which reduce restitution slope, play a protective role against arrhythmias. Accordingly, computational experiments can effectively support the design of new experiments with hiPSC-CMs.

Conclusions
Computational simulation of cardiomyocytes' electrophysiological properties may help to explain the shape of the electrical restitution curves at various repolarization levels for atrial and ventricular cardiomyocytes and to predict the kinetics of intracellular calcium and contractile force. The proposed model allowed us to detect those regions characterized by elevated slopes, capable to confer arrhythmogenic properties of cardiomyocytes. Computational results are useful to interpret experimental results with hiPSC-CMs on the lab-on-chip platform and to propose the new design of the experiments for personalized studies of heart disease. The aforementioned research can allow monitoring the changes in AP, calcium current properties of hiPSC-CMs and the development of arrhythmias in response to application of various drugs or different stimulation modes. The method of dynamic analogies used in the work allows us to generalize this approach for studying the electrical properties of cardiomocytes and for investigation of the electromechanical model of nonequilibrium processes in the myocardial tissues. Future directions of research related to the simulation of cardiomyocyte's activity should be focused on electro-mechanical coupling, electrical and mechanical stimulation, and spontaneous beating