Numerical Simulation of Excitation-Contraction in Isolated Cardiomyocytes

We present and study a mathematical model to simulate the calcium ﬂux and contractile activity of the cardiomyocyte, resorting on the Finite Element Method. This Model of cardiac cells provide the possibility to understand the biochemistry and biomechanics of cardiac cells and cardiovascular diseases. Heart Failure is considered the ultimate cardiac disease, a condition with no effective cure which is highly related with the function of cardiomyocytes and consequently with the concentration of calcium, affecting the contractile activity of the cardiomyocyte. In order to test our model’s performance, several tests were applied varying the local active cellular tension driven by the intracellular calcium concentration and the localization of the main calcium inﬂux. The results are expressed in the graphics of calcium concentration over time, maximum cardiomyocyte contraction and the gradient of calcium diffusion. Our results show that the behaviour of our models is faithful to what is known to be true with cell’s physiology and pathological conditions.


Introduction
The recent evolution in biotechnology and informatics has promoted the production of large amounts of information and has accelerated the processes of investigation in biological systems and cellular-based studies.The use of mathematical and geometric models enables the simulation of biological complex systems, and therefore, the formulation of hypotheses and realisation of experiments that would normally require the use of live cells and its lifetime conditions.
Models of cardiac cells are capable of facilitating insights into the mechanisms underlying cardiac dynamics and its pathophysiology.They provide the possibility to understand the biochemistry and biomechanics of cardiac cells, since the cardiovascular diseases have reached an enormous burden on contemporary society.In fact, Heart Failure, the ultimate cardiac disease which has no effective cure, has been increasing in mortality and morbidity [1], and its evolution always has a reserved prognosis.From this very negative trend, the need arises to create and test a new tool that could be an added value in the study of the biomechanical mechanisms of the cardiomyocyte: numerical modelling.[2] In this work, the numerical modelling uses equations, discretized by the Finite Element Method, which model the known cellular dynamics, as well as the possible dysfunctions.Thus, it becomes possible to create and to experiment with other hypotheses, with unlimited resources, and without raising ethical questions.
The human heart is an essential organ with a highly organized tissue constituted of ventricular and atrial cardiomyocytes, pacemaker cells, Purkinje cells, vasculature, and connective tissue.Cardiomyocytes are mononuclear cells and form the striated muscle of the heart, with the shape of a tridimentional net.This tridimentional structure allows the heart to contract and relax in the most efficient way.The ventricular cardiomyocytes are columnar shaped cells of 20nm in diameter and 60 − 140nm in length, while the atrial cardiomyocytes are ellipsoidal shaped cells of 5nm in diameter and 10 − 20nm in length.[3] Figure 1.Morpholohy of Cardiac Muscle [4] The biggest part of the whole cell volume (75%) is made up of myofibrils and mitochondria.Myofibril is the bundle that forms the contractile elements within the cell and it is composed of actin thin filament, myosin thick filament and titin, which stabilizes myosin at the Z-line.The fundamental unit within the cardiomyocyte that allows its contraction is the sarcomere (aggregation of myofibrils), which has a length of 1.8µm in the systole and 2.2µm in the diastole.
Other than myofibrils, this contractile structure contains tropomyosin (troponin complex).Myosin has a filamentous tail and a globular head region that contains the site for actin binding.Actin has F-actin, the backbone of the thin filament, and G-actin, which works as a stabilizing protein.
The Sarcoplasmic Reticulum (SR), the nucleus, and the cytosol occupy the other 25% of the cell.Surrounding these components, there is a specialized structure called the Sarcolemma, which is a coalescence of the plasma membrane and the basement membrane, composed of a lipid bilayer.This wall contains hydrophilic heads and hydrophobic tails, allowing the regulation of the interactions with the intracellular and extracellular environment.In order to maximize the in and out interactions, there are transverse tubular systems (T-tubules) in the sarcolemma.
The T-tubules are 150nm to 300nm wide deep invaginations of the sarcolemma into the cardiomyocyte, and they form an extended barrier between the intracellular and extracellular space, turning them into important structural components in the Excitation-Contraction Coupling (ECC) system.

NUMERICAL SIMULATION OF EXCITATION-CONTRACTION IN ISOLATED CARDIOMYOCYTES
All of these cell organelle interactions are included in the ECC process and are fuelled by the flow of three important ions: K + , Ca 2+ and N a 2+ .

Structural details
The Ryanodine Receptor (RyRs) are ion channels that open in the presence of elevated Ca 2+ and release calcium from the junctional sarcoplasmic reticulum.The junctional sarcoplasmic reticulum makes such close contact with the T-tubule membrane that RyRs channels on the sarcoplasmic reticulum are impending apposed (≈ 15nm) to L-type Ca 2+ channels on the T-tubule (≈ 14 in peripheral couplings to ≈< 100 in intracellular sites), this way forming the cardiac dyad that is fundamental to the processes of Ca 2+ flow.from the junctional sarcoplasmic reticulum.When viewed on detail, T-tubules are seen to radiate throughout the cardiac cell.This close association between the T-tubule and sarcoplasmic reticulum ensures the synchronous rise of [Ca 2+ ], during systole, therefore, the importance of the whole architecture of the cell.In addition to these basic requirements for ECC process, the cardiac dyad may also be considered as containing additional structures that may contribute to or modulate Ca 2+ release from the sarcoplasmic reticulum during systole.The most extensively studied is N a 2+ /Ca 2+ eXchanger (NCX), that has been argued via its reverse-mode action to contribute to Ca 2+ influx early during the action potential.The second protein of interest is SERCA2a: Sarco(Endo)plasmic Reticulum Calcium ATPase.Available data on SERCA2a distribution in the heart is surprisingly sparse; however, confocal studies seem to suggest that a significant proportion of SERCA is localized to the z-line as well as enveloping the myofilaments.[5] This cycle of calcium triggering, release and uptake, takes place in the approximately 20, 000 diad spaces (also called calcium release sites) that are distributed throughout the myocyte.[6] Figure 2. Close look to cardiac diad [7] 3. From Action Potential to Contraction ECC characterizes the whole processes relating to electrical excitation through force generation and the power of contraction in the heart.It occurs at multiple levels from the whole heart, to each single cardiomyocyte and down to the sarcomere.[6] The propagation of the action potential through the entire cardiac muscle requires the depolarization of a cardiomyocyte.This means that the cell will generate an action potential while arising from a steady state to an excited state, and then pass it to the neighbour cell.In order to induce the process, an electric charge will be propagated by the three cations already referred.These three ions are small enough to be widespread through the cells membranes.The contraction of individual myocytes is coordinated by a propagation wave of electrical depolarization of the sarcolemma.
The propagation from cell to cell of the action potential corresponds to the beginning of the action potential cycle (phase 0) and requires a difference between cell voltages.The bigger the difference, the bigger the flow between the cells.When the stimulus reaches a steady state cell, it becomes excited and some of the N a 2+ channels open.The N a 2+ moves in the cell, becoming less negative.When the cell voltage reaches −70mV , every N a 2+ channels open and more N a 2+ ions get in the cell.The cell becomes positive.The N a 2+ channels close, making the cell more negative again.607

Chemical aspects
In the beginning of the contraction stimulus, the K + channels open and K + ions move out of the cell.This K + extraction, decreases the cell potential from 20mV to 0mV (phase 1).Ion Ca 2+ is fundamental for evoking the ECC complex.Upon the wave of depolarization, when electrical action potential reaches T-tubules, the wave of depolarization induces Ca 2+ influx into the cardiomyocyte through the voltage-sensitive L-type Ca 2+ channel of the T-tubules.This rapid but small Ca 2+ influx causes activation of large amounts of Ca 2+ release from the RyR2 on the sarcoplasmic reticulum.Finally, cytosolic Ca 2+ level changes from 100nmol/L to 10µmol/L in concentration.The 10µmol/L of Ca 2+ binds to the troponin C, leading to the actin-myosin interaction resulting in initiating crossbridge formation and, therefore, contraction movement (phase 2).
Using Adenosine Tri-Phosphate (ATP), the G-actin interacts with the myosin globular head leading to the crossbridge formation and sarcomere shorting (contraction).Tropomyosin lies on the side of actin for rigidity to thin filament.The troponin complex, also present in the thin filament, is composed of troponin C, I and T.These proteins are responsabile for the regulation of the crossbridge formation.
Active relaxation (phase 3) of the cardiomyocyte is dependent on the function of the SERCA2a.For each 1mol of ATP hydrolysed, 2mol of Ca 2+ is transported back into the sarcoplasmic reticulum.Phospholamban (PLB) regulates the function of SERCA2a.Additionally, the NCX on the plasma membrane removes Ca 2 from cytosol, until the cell reaches a voltage of −90mV , changing Ca 2+ for N a 2+ .In summary, the process of ECC links the electric excitation of the surface membrane, i. e. action potential, to contraction, with the concentration of ion Ca 2+ being the major factor for the magnitude and duration of the contraction.[9] In the end, the purpose of cardiomyocytes is to generate a controlled contraction in response to repetitive electrical depolarization signals -a function that is fully controlled by Ca 2+ .
The importance of ion Ca 2+ will be shown throughout this work, as the ion's flow will be used as an effective tool to understand the complex dynamics of the processes involved in excitation-contraction coupling and in its computational simulation.

Mathematical model aspects
We consider for a seek of simplicity the quite well established two variable model of Goldbeter et al. [11] as the core for self-sustained calcium oscillations (phases 2 and 3).In this model, the calcium-induced calcium release mechanism is described by nonlinear kinetics between cytosolic c = c(t, x, y) and sarcoplasmic s = s(t, x, y) calcium concentrations at instant t and position (x, y), with rates of variation given by the differential system: 608 The intracellular propagation of free calcium is taken into account in the model by considering the anisotropic diffusion of calcium observed experimentally at [13] with In this model, v 1 β is the constant flow of calcium into the cytosol, controlled by the sarcoplasmic reticulum.The concentration of free cytosolic calcium c oscillates between low levels, during which c is pumped by the sarcoplasmic reticulum (v 2 flux), and high levels, characterised by the autocatalytic release of calcium v 3 from the sarcoplasmic reticulum.The flux k s s is a basal leak of calcium into the sarcoplasmic reticulum, while the input flux v 0 and efflux k c c refer to the Ca 2+ fluxes into and out of the cell, respectively, finally D is the diffusion tensor with components D ij .
Equations ( 1) and ( 2) describe an entire family of models dependent on the particular choice of the exchange function F (c, s) that describes the dynamics of Ca 2+ exchange between the cytosol and the internal Ca 2+ pool.It should, therefore, be kept in mind that, although we present results here for a particular choice of F , similar results can be obtained from exchange functions with similar qualitative properties.
For the contraction movement (phase 2) we assume that the total stress σ = σ ij within the cardiac cell is generated by the contraction of actin-myosine fibers.Since the sarcomeres are oriented along the cell principal axis (x axis), we consider an isotropic linear elastic motion governed by: where γ (c) T u is the local active cellular tension driven by the intracellular calcium concentration c = c(t, x, y).
We assume that the analytical expression of γ is given by a Hill function of the form [14]:
We remark that at equation (10) we have used the Hooke's law for plane stress, with the displacement vector (u x , u y ) for the relationship between the two-second-order tensors, the strain tensor ε and the stress tensor σ: σ = λtr(ε(u x , u y )II + µε(u x , u y ), (11) where II denotes the 2-identity matrix, and tr(•) represents the trace of a matrix.The Lamé positive constants λ and µ are given by (cf [16]) with the Young's modulus E > 0 and the Poisson ratio ν ∈ ]0, 1/2[.
The time integration will be done by using a backward Euler method which is written as: (15) for the instant τ > 0 and a given guess ,where h is called the step size.We remark that at each instant we obtain a new configuration for the cell due to the displacement field ⃗ u τ +1 = (u τ +1 x , u τ +1 y ) coming from the concentration c τ +1 .We refer [15] for theoretical questions about existence and unicity of the numerical solutions over adequate finite element spaces.
By reasons of numerical stability the equation ( 13) must be penalized with an artificial diffusion therm using a small parameter ξ like that: The numerical scheme can be defined by:    Given a given guess , τ > 0 by solving eqs( 16), (14) and (15) resp.(17) We consider the cycle ends when the quantity of calcium is low in the cell: ∫ Ω c τ +1 dΩ < ϵ , for a small ϵ.

Numerical results
At this section we will apply scheme (17) with the parameters from Table 1 using the partial differential equation solver, using finite element methods, FreeFem++ [17] .
Table 1.Cell parameters values [11]  For all cases we simulate only one cycle, supposing a conveyable conditions of ATP and oxygen.To better explain the results, we will present the graphic of calcium concentration over time, maximum cardiomyocyte contraction and the gradient of calcium diffusion for each case.
For all exhibited graphics we have used the time τ as unity for the x-axis and in the y-axis we will consider calcium concentration µM for the blue graphic and percent (%) of shortening for the green one.

Case 1
In this case, as described before, T-tubule is placed on the left side of the cell, and the right side is free to contract.The cardiomyocyte's contraction happens when the calcium concentration is greater than µM, approximately, and we remark with Figure 5 that the contraction peak time is approximately in one-third of the cycle as it was expected.The maximum cardiomyocyte contraction is achieved at τ = 41 ms (millisecond), as we show in Figure 6.The initial cell length is 10 µm, and after contraction its length is about 8.9 µm, which results in a reduction of 11% (% of shortening), s shown at Figure 6.The simulation of the T-tubule effect is exhibited in Figure 7 with the gradient of the calcium concentration c τ , τ = 41 ms.

Case 2
In this case 2, T-tubule is also placed on the left boundary of the cell, but both left and right sides are free to contract.The cardiomyocyte's contraction happens when the calcium concentration is greater than 0.45 µM, approximately, and we can remark with Figure 8 that the contraction peak time is approximately one-third of the cycle as it was theoretically expected.The maximum cardiomyocyte contraction is achieved at τ = 41 ms, as we show in Figure9.At this case the cell length after contraction is about 9.4 µm, which results in a reduction of 6%.

Case 3
In this case 3, T-tubule is placed on the right side of the cell, and both left and right sides are free to contract.The cardiomyocyte's contraction happens when the calcium concentration is greater than 0.45 µM, approximately, and we can remark with Figure 11 that the contraction peak time is in approximately one-third of the cycle as it was expected.The maximum cardiomyocyte contraction is achieved at τ = 41 ms, as we show in Figure 12.As above the initial cell length is 10 µm, and after contraction its length is about 9.3 µm, which results in a reduction of 7%.The simulation of the T-tubule effect is exhibited in Figure 13 with the gradient of the calcium concentration c τ , τ = 41 ms.

Case 4
In this case 4, T-tubule is placed, over a reduced length zone, on the top and bottom boundaries of the model (Figure 4 d) ) and both left and right sides are free to contract.The cardiomyocyte's contraction happens when the calcium concentration is greater than 0.36 µM, approximately, and we can remark with Figure 14 that the contraction peak time appears in the middle of the cycle.
The maximum cardiomyocyte contraction is achieved at τ = 82 ms, as we show in Figure 15.After contraction the cell length is about 9.8 µm, which results in a reduction of 2% .
The simulation of the T-tubule effect is exhibited in Figure 16 with the gradient of the calcium concentration c τ , τ = 82 ms.
Finally we resume at Table 2 the main results from the study cases.

Conclusion
We have investigated an Finite Element Method approach for simplified coupled mathematical model of excitationcontraction in isolated cardiomyocytes.We started with the variational formulation model and proceed its numerical implementation based on software FreeFem++.
The results from the study cases highlight to the T-tubules distribution importance along the cardiomyocyte.We remark that for all cases the cell deformation is higher near the T-tubule boundary, as the calcium concentration.Table 2 points to the importance of a large zone for calcium influx and a existence of sides free to contract, this simulate the contact between two cells.
For case 1, the result was a concentration of 0.45µM and a % of shortening around 10% corresponds to the results obtained by DelMonte et al [18] and Høydal et al [19].These results are the ones expected for the physiological function of a real cardiomyocyte.This means that this simulation is reliable for these two parameters, concentration and % of shortening.
In case 2, the result was a concentration of 0.45µM and a % of shortening around 5% corresponds to the results obtained by Kamgoue et al (3%) [20].These results can also be considered normal or indicate a situation of hypocalcemia (3%).
The results of case 3 were presented with a concentration of 0.45µM and a % of shortening from 4% to 10% , which corresponds to the results obtained both by Kamgoue et al [20] and Garcia-Canadilla et al (9%) [21].
At last, for case 4, with the special location of T-tubules (on the top and bottom boundarys), we reach a concentration of cytoplasmatic calcium of 0.36µM and a % of shortening of 2%, which means that for this kind of model distribution (T-tubule on top and bottom), the local active cellular tension must to be higher for the same % of shortening.
The results obtained with our model are faithful to what is known to be true with cell's physiology and pathological conditions (hypocalcemia).
This study proposes a simplification of the mathematical context of the excitation-contraction in isolated cardiomyocytes problem and brings a contribution to a very complex subjet.
A relevant study is to introduce in this mathematical model multiple parametric, to enable understanding its influence, namely the ones that are considered to be important to studied cardiac diseases.

Figure 3 .
Figure 3. Action Potential in Cardiomyocytes

× 10 4
PaWith the purpose of fulfilling a wide range of model cases, we created four different models of cardiomyocytes, varying the possible zones to place the mathematical boundary conditions for the simulation of the T-tubules, as well as the freedom of each side of the cell, in order to obtain different kinds of contraction movements.By symmetric reasons, we describe the four case as the following: 1. Case 1: T-tubule is placed on the left side of the cell, and the right side is free to contract, Figure 4 a).2. Case 2: T-tubules is placed on the left side and both sides are free for the contraction movement, Figure 4 b).3. Case 3: T-tubule is placed in the right side and both sides are free for the contraction movement, Figure 4 c).4. Case 4: Both sides are free for the contraction movement, but there are T-tubules from the upper boundary and from the bottom boundary , Figure 4 d).

Figure 4 .
Figure 4. Outline of cardiomyocyte with different localizations for Γ T

Figure 8 .
Figure 8. Cytosol concentration in the cell for one cycle, in Case 2 32

Figure 11 .
Figure 11.Cytosol concentration in the cell for one cycle, in Case 3 35

Table 2 .
Main results comparison for study cases