Modeling the cardiac electromechanical function: A mathematical journey

By Alfio Quarteroni, Luca Dedè, and Francesco Regazzoni

Abstract

In this paper we introduce the electromechanical mathematical model of the human heart. After deriving it from physical first principles, we discuss its mathematical properties and the way numerical methods can be set up to obtain numerical approximations of the (otherwise unachievable) mathematical solutions. The major challenges that we need to face—e.g., possible lack of initial and boundary data, the trade off between increasing the accuracy of the numerical model and its computational complexity—are addressed. Numerical tests here presented have a twofold aim: to show that numerical solutions match the expected theoretical rate of convergence, and that our model can provide a preliminary valuable tool to face problems of clinical relevance.

1. Introduction

Historical tips

Since the dawn of history, the heart and its functioning have fascinated the best minds on our planet, through a long journey marked by intuitions and dazzles. Three centuries bc, Aristotle claimed that the blood vessels transmitted animal heat from the heart to the periphery. Shortly thereafter, Praxagoras of Kos sensed the different role of arteries and veins and assumed that the arteries transported air, and veins blood. It was Galen in the second century ad to recognize that arteries also carry blood.

Thirteen centuries later, the great Leonardo Da Vinci distinguished for the first time ventricles from atria (auriculae), describing the moti del core as a result of the contraction of the atria that coincides with the ventricular diastole, while, conversely, the ventricular systole occurs along with the expansion of the atria. In the seventeenth century Harvey correctly described the blood circulation for the first time. Thanks to the beating heart vivisection of various animals, he recognized that the heart valves in the veins were designed to activate only if the blood flowed to the heart and therefore concluded that “the blood has a movement, and it is circular”. In the following century, the great mathematicians Leonhard Euler and Daniel Bernoulli made decisive steps forward. In 1730 Bernoulli, professor of mathematics and anatomy at the University of Basel, formulated his famous vis viva equation between pressure, density, and velocity of blood. In 1775 Euler, in a work entitled Principia pro motu sanguinis per arteria determinando Reference 30, wrote a series of differential equations—still fundamental and used in areas as diverse as the gas motion in pipes or aircraft aerodynamics design—to describe the evolution of blood flow and pressure in an idealized cylindrical blood vessel. These equations have then been generalized and today still set the ground for a one-dimensional model of blood flow in the complex network of arteries and veins of our circulatory system.

An integrated heart model

The ambition of today’s mathematicians is to travel the furrow traced by Euler for creating a virtual, immaterial heart, made only of equations, able to reproduce all the vital cardiac processes: the propagation of the electric potential, the contraction and relaxation of the myocardium (the heart muscle), the blood dynamics in atria and ventricles, the coronaries perfusion of the myocardium, the valve dynamics, and all their mutual interactions. This extraordinary complexity, which mysteriously determines the harmonious synchronization producing the heartbeat, can be turned into a mathematical model, thanks to the first principles of physics, that is Newton’s laws, thermodynamics, chemical kinetics, Navier–Stokes for blood dynamics and Maxwell’s for the propagation of the electric field, etc. However, this is not enough, as the corresponding equations need to be supplemented with constitutive equations that characterize the specific nature of the cardiac tissue at hand. One such relationship relates stresses and strains for the cardiac muscles and accounts for the material properties of the myocardium. Other equations are necessary to characterize the way cardiomyocytes (the elementary components of the cardiac muscle) contract and relax once stimulated by electrical signal propagating through the myocardium. The result is a system of coupled nonlinear differential equations, ordinary differential equations (ODEs), and partial differential equations (PDEs): once properly solved—and this will inevitably call into play computers, as a matter of fact supercomputers!—this system ought be suitable to simulate the human heart function. The solution of this integrated heart model (IHM) represents a scientifically arduous challenge: the dream is however to provide researchers in applied mathematics, bioengineering, and life sciences, as well as physiologists, cardiologists, and cardiac surgeons, with an effective and powerful tool aimed to the qualitative and quantitative study of the heart, both in physiological or pathological conditions, and at making advances in the diagnosis and treatment of heart diseases.

In this paper, we provide an in-depth and rigorous mathematical presentation of the most crucial aspect of the entire IHM: the so-called electromechanical (EM) mathematical model. We will first derive it from basic principles. In particular, we will highlight the way the activation force triggered by the transmembrane electric potential is generated at the cardiomyocytes level and then affects the whole muscle deformation. We will then recall the main mathematical results inherent to the well-posedness of the electromechanical model components.

The role of data

A mathematical model, if properly designed, enjoys the property of universality, that is its validity is independent of the specific context to which it is applied. Of course, what makes the mathematical model most valuable is its ability to adapt to any specific context: in our case, virtually to every possible individual, thanks to the data that characterize it. Data are needed to feed models. For an IHM, data concern the geometrical shape of a specific heart, initial and boundary conditions that must supplement the differential equations, parameters that characterize material properties, just to mention a few. Data is a crucial issue. Shapes can be retrieved from medical images (such as computer tomography or magnetic resonance imaging) after a good deal of geometric preprocessing. Initial and boundary conditions are seldom available, as they would require invasive (and unnecessary) clinical examinations on patients. Also material parameters are not easy to obtain for specific patients due to the difficulty of producing ad hoc measurements. To compensate for the insufficiency (or lack) of data, mathematicians have historically came up with brilliant ideas such as parameters identification, variational or Kalman filter-based techniques for data assimilation. Moreover, whereas the mathematical model is deterministic in itself, data could be uncertain. It is therefore of paramount importance to quantify how this uncertainty propagates from the input data to the outputs. Dealing with uncertainty in data has required ad hoc, often sophisticated, mathematical tools as well.

Road to solution

Once the differential system is completed by the missing data, closed-form solutions seldom (almost never) exist. The first step is an a priori mathematical analysis to understand in which regimes a solution exists, whether it is unique and if it continuously depends on data (this is the well-posedness in the sense of Hadamard).

Once our problem is set upon rigorous mathematical ground, the next step consists of approximating it by reducing it to a finite dimension, say $N$ (typically very large, of an order of millions or more). A mathematically rigorous strategy consists in projecting the variational formulation (that is the weak formulation) of the differential system onto finite-dimensional subspaces of the Sobolev spaces (typically $H^1$, see, e.g., Reference 3) to which the weak solution of the original system belongs. To make sure that we have suitably operated, our new finite-dimensional (numerical) problem must enjoy the properties of stability and convergence. Stability is the continuous dependence (uniformly, with respect to $N$), in suitable norms, of the solution to the problem’s data. Convergence is the property of the numerical solution to tend to the exact (unknown) one once $N$ tends to infinity. Error estimates represent a further peculiar feature of numerical problems. Besides being convergent, the numerical solution should converge to the exact one in suitable norms with a specific rate of convergence with respect to $1/N$. In practice, for models of this kind of complexity, a tradeoff between accuracy and computational effort has to be pursued. Needless to say, as $N$ is unavoidably extremely large (often reaching hundreds of millions), supercomputers are necessary at this stage. Indeed, the simulation of a single heartbeat (about one second of physical time) may require several hours of computation on a big, modern supercomputer.

But what is all that for? From a mathematician’s perspective, mathematical and numerical models of the heart function represent an inexhaustible source of mathematical challenges. From an applied mathematician’s perspective, the excitement is further enhanced by the ambition to provide doctors with a “tool” that may help them with diagnosis and treatment.

In this paper we will illustrate the paradigm (data, mathematical model; numerical approximation, computer simulation), and we will show in a few examples the way we can nowadays provide doctors with solutions of clinical relevance. In this regard cardiac electromechanics provides an excellent testing ground, more and more as the state of the art in mathematical modeling, numerical methods, and computational platforms advances. As a matter of fact, starting from the visionary works of C. S. Peskin Reference 62Reference 63, this topic has inspired an impressive body of research from many research groups around the world (see, e.g., Reference 10Reference 11Reference 17Reference 21Reference 25Reference 26Reference 32Reference 37Reference 39Reference 57Reference 65Reference 83Reference 85Reference 86Reference 90Reference 94).

Paper outline

The outline of the paper is as follows. In Part 1 we will give a short introduction to the elementary concepts of cardiac anatomy and physiology. Then we will introduce the mathematical model of cardiac electrophysiology and active force generation. We will then focus on the modeling of cardiac mechanics and will conclude with the multiscale model of cardiac electromechanics. Part 2 will be devoted to the finite-dimensional approximation in space (by the finite element method) and in time. We will then comment on the use of efficient solvers for the associated nonlinear algebraic systems of very large size. Part 3 will then illustrate the numerical solutions that we can achieve, as well as some examples that are relevant to clinical applications. Conclusions and possible generalizations follow.

Part 1. From the physical problem to the mathematical model

2. Cardiac anatomy and physiology at a glance

The human heart is a sophisticated machine, whose functional role is pumping blood throughout the body to its cells, supplying the organs with oxygen, and removing the metabolic waste. The contraction of the heart muscle is due to that of the cardiomyocytes, the cardiac cells, once activated by electrical processes. As a matter of fact, at each heartbeat, an electrical potential propagates through the heart tissue until it reaches each cardiomyocyte of the atria and ventricles (the heart chambers). This synchronized propagation activates an orderly contraction of the heart chambers and generates the heart rhythm. More specifically, each heartbeat is triggered by an electrical signal, originating from the sino-atrial node, the heart’s natural pacemaker consisting of a group of self-excitable cells that is located in the upper part of the right atrium (see Figure 1). The signal propagates from one cell to another through the two atria and reaches the atrioventricular node, located between the atria and the ventricles. The atrioventricular node acts as a filter for signal propagation to ensure that the contraction of the ventricles begins once the blood has passed from the atria into the same ventricles. The electrical signal then propagates from the atrioventricular node through the Purkinje fibers, a high-conductibility network of fibers, reaching the myocardial cells.

In order to derive a mathematical description of these macroscopic processes, it is necessary to start from a description of the microscopic mechanisms that generate them, following a micro-to-macro approach. This process starts from the elementary components of the cardiac muscle tissue, cardiomyocytes (i.e., cardiac muscle cells). These are in fact excitable: when electrically stimulated, the electrochemical balance of the cell membrane changes, giving rise to a sequence of biochemical processes that determine a significant variation in cell potential: rapid depolarization followed by a repolarization. This phenomenon, known as action potential, is due to the opening and closing of ion channels, located in the cell membrane. The latter becomes permeable to different ions (calcium, potassium, magnesium) thanks to the transmembrane potential, that is the difference in voltage between the internal and external part of the cell. The ionic fluxes determine a variation of the transmembrane potential and have a feedback effect on the voltage difference itself. Among the various ionic species involved in the dynamics of the action potential, calcium ions play an important role. Calcium represents the trigger for muscle contraction: calcium ions induce a complex chain of reactions generating an active force inside the cardiomyocytes. Finally, thanks to a process of transmission through different space-time scales, the active force at the microscopic scales (that of cardiomyocytes) generates a resultant force at the organ level, thus giving rise to the contraction of ventricles and atria. The combined effect of active and passive force (i.e., the reaction of the myocardium to mechanical stress), as well as the coordinated action of the atria and ventricles, governs the blood fluid dynamics in the four chambers as well as the valve dynamics. The building blocks of this complex network of interactions are summarized in Figure 2. In Sections 3, 4, 5, and 6, starting from the microscale and progressively climbing the hierarchy of scales up to the macroscale, we present mathematical models describing the above mentioned processes.

3. Modeling cardiac electrophysiology

As mentioned, the driver of the cardiac function is electrophysiology, i.e., the result of chemical and electrical processes taking place at different spatial scales, from subcellular to the the whole organ scale. A mathematical description of these processes is based on the translation into mathematical terms of the principles ruling the electro-chemical activity of ions species at the finest scale; then, by progressively climbing up the hierarchy of spatial scales, a set of equations describing the tissue-level electrophysiological activity is derived.

Let $\Omega _{0} \subset \mathbb{R}^3$ be an open connected set, denoting the region of space occupied by the cardiac tissue. This region can be split into two subregions, namely the inner space of cells (intracellular space) and the region located outside the cell membranes (extracellular space). The intracellular and extracellular spaces are separated by the cell membranes, which are scattered by ion channels. An ion channel can be understood as an opening in the cell membrane that selectively allows the transit of specific ion species. The movement of electrically charged ions generates an electric current. Starting from Maxwell’s laws of electromagnetism and Einstein’s theory on Brownian motion Reference 29, one can derive a system of PDEs describing the ion dynamics. An analytical solution of this system can be obtained under the (realistic) assumption that the ion dynamics in the direction tangential to the cell membrane is negligible compared to the transmembrane one. By further neglecting the space variations across the membrane and focusing exclusively on space-averaged equations, this allows us to link the electric current crossing the membrane (denoted by $I_{\mathrm{ion}}$) with the concentration of ionic species and with the transmembrane potential, defined as ${v} = {u_{i}} - {u_{e}}$, where ${u_{i}}$ and ${u_{e}}$ denote the electric potential in the intracellular and extracellular spaces, respectively.

The membrane of cardiomyocytes is selectively permeable to multiple ionic species, such as sodium, potassium, and calcium, whose fluxes are regulated by the opening and closing of the ion channels. To track the opening and closing of these channels, as well as the resulting electric current, several models have been proposed in the literature. These models are written in the general form

$$\begin{equation} \frac{d \mathbf{z}_{\mathrm{ion}}}{dt} = \boldsymbol{\Phi }_{\mathrm{ion}}({v}, \mathbf{z}_{\mathrm{ion}}), \qquad I_{\mathrm{ion}} = \mathscr{I}_{\mathrm{ion}}({v}, \mathbf{z}_{\mathrm{ion}}), \cssId{texmlid1}{\tag{3.1}} \end{equation}$$

where $\mathscr{I}_{\mathrm{ion}}$ and $\boldsymbol{\Phi }_{\mathrm{ion}}$ are suitably defined functions and $\mathbf{z}_{\mathrm{ion}}(t) \in \mathbb{R}^{N_{\mathrm{ion}}}$ is a vector collecting the so-called ionic variables, describing ion channels, concentrations, or simply phenomenological variables. For the latter, a typical example concerns recovery variables, namely variables that account for the tissue refractoriness, without a one-to-one relationship with measurable physical quantities. Ideally, system Equation 3.1 applies at any cell location. Notable examples are provided by the Aliev–Panfilov model Reference 4, which envisages a single ionic variable (i.e., $N_{\mathrm{ion}} = 1$) that accounts for tissue refractoriness; the Ten Tusscher–Panfilov model Reference 89, which includes 18 variables describing ionic concentrations and the probability of opening of specific ion channels; the O’Hara–Rudy model Reference 59, featuring $N_{\mathrm{ion}} = 41$ variables, which was built and calibrated upon data coming from more than 100 undiseased human hearts.

To complete a model describing the electrical activity across the cell membrane, we need to relate the total ionic current $I_{\mathrm{ion}}$ to the dynamics of the transmembrane potential ${v}$. With this aim, the cell membrane is modeled as a capacitor with capacity ${C_{m}}$: when a current $I$ crosses the cell membrane for an infinitesimal time interval $dt$, the transmembrane potential difference decreases by $I / {C_{m}} \, dt$ (by convention, current and potential difference are taken with opposite signs).

Furthermore, equation Equation 3.1 describes the electrical activity of a single cell. However, in the cardiac tissue, the nearly three billion cells are not electrically isolated. Indeed, they are connected by the so-called gap junctions, which link the intracellular and extracellular spaces of adjacent cells and which allow the electrical signal to propagate from one cell to another. The following equation, known as the monodomain model Reference 24, accounts for these effects by means of a diffusive term:

$$\begin{equation} \chi _{m} {C_{m}} \partial _t {v} - \nabla \cdot \left( \mathbf{D} \, \nabla {v} \right) + \chi _{m} I_{\mathrm{ion}} = \mathcal{I}_{\mathrm{app}} \qquad \text{in } \Omega _{0} \times (0,T]. \cssId{texmlid2}{\tag{3.2}} \end{equation}$$

In the above equation $\mathcal{I}_{\mathrm{app}}$ represents an externally applied current, while $\chi _{m}$ denotes the area-to-volume ratio (i.e., the amount of membrane area per volume of tissue), and it allows us to relate the microscopic cell dynamics to quantities related to a macroscopical characterization of the tissue. The tensor $\mathbf{D}$, known as the conductivity tensor, characterizes the electrical properties of the tissue. A feature of cardiac tissue is that the electrical signal propagates with a larger speed in the fibers’ direction (i.e., $\mathbf{f}_{0}$) than in the orthogonal directions. To account for such an anisotropic property, $\mathbf{D}$ is typically defined as

$$\begin{equation} \begin{split} \mathbf{D} &= \sigma _{f} \mathbf{f}_{0} \otimes \mathbf{f}_{0} + \sigma _{s} \mathbf{s}_{0} \otimes \mathbf{s}_{0} + \sigma _{n} \mathbf{n}_{0} \otimes \mathbf{n}_{0} , \end{split} \cssId{texmlid6}{\tag{3.3}} \end{equation}$$

where $\sigma _{f}$, $\sigma _{s}$, and $\sigma _{n}$ represent the conductivity coefficients in the fibers, sheets, and cross-fibers directions. The triplet $(\mathbf{f}_{0},\mathbf{s}_{0}, \mathbf{n}_{0})$ not only characterizes the electrical behavior of the tissue, but also its mechanical response, as we will describe in Section 5.

The applied current $\mathcal{I}_{\mathrm{app}}$ is typically prescribed as a function with support restricted to a short time interval (usually a few milliseconds) at the beginning of each heartbeat and to a few spots located close to the surface of the myocardium, representing the end points of the Purkinje fibers. The typical solution of Equation 3.2 is in the form of a traveling wave for the variable ${v}$, originating from the locations where $\mathcal{I}_{\mathrm{app}}$ is applied. When a point of the domain is reached by the wavefront, ${v}$ quickly raises (depolarization), then features a plateau, and finally returns to its resting position (repolarization); this yields the so called action potential. Due to the quick depolarization, the potential wave is characterized by a very steep front.

The monodomain model can be derived from the bidomain model, in which two separate variables are associated with the electric potential of the intracellular and extracellular space. In its turn, the bidomain model can be rigorously derived by assuming the existence of the extracellular and intracellular spaces as two simply connected subdomains and by employing a homogenization technique Reference 24.

A general comment is in order. Even when it is derived from physical principles, the mathematical meaningfulness of a model is not guaranteed a priori yet. Specifically, one should check that the problem is well-posed from a mathematical standpoint, that is, a solution exists, it is unique, and it depends continuously on data (the importance of the latter assumption is related to the uncertainty that unavoidably affects the measurements of these data). Moreover, the existence of solutions is an essential requirement to address the numerical approximation of a model. As a matter of fact, even if approximate solutions exist, they would be meaningless in the absence of an exact solution to which they can converge. In this regard, we report that a well-posedness result for the bidomain model coupled with the FitzHugh–Nagumo model Reference 33 is presented in Reference 34. Another well-posedness result, that makes use of a fixed-point argument, is provided in Reference 93.

4. Modeling active force generation

Due to the action potential dynamics, in the first stages of each heartbeat the calcium concentration inside cardiomyocytes (denoted by $[\mathrm{Ca}^{2+}]_{\mathrm{i}}$, where i stands for intracellular) quickly raises by nearly one order of magnitude before returning to its resting value in nearly a half second. In fact calcium acts as a cellular messenger by triggering the contraction of cardiomyocytes. At the microscopic scale, cardiomyocytes are organized in sarcomeres, cylinder shaped contractile elements of the size of nearly 2 μm. Sarcomeres are composed of myofilaments, thin filaments made of the proteins troponin, tropomyosin, and actin, and thick filaments made of the protein myosin. The active force is generated by the interaction between actin and myosin, molecular motors capable of transforming the chemical energy into mechanical work Reference 47. The long-standing hypothesis that muscle contraction was led by folding of elongated protein filaments was challenged by the discovery that the filaments’ length remains constant during contraction and that it is instead the mutual sliding between the two families of filaments (thin and thick) that makes the muscle contract Reference 45Reference 46. This discovery, known as sliding filaments theory, was made independently by two research teams: on one side, the British biologist H. Huxley and biophysicist J. Hanson, working at MIT; on the other, the British physiologist A. F. Huxley (Nobel prize winner in 1963 for his work on the action potential) and the German physician R. Niedergerke, working at the University of Cambridge. The two teams decided to publish their work in two consecutive articles in the same issue of Nature Reference 45Reference 46.

Many models describing these microscopic mechanisms have been proposed in the literature. In general terms, these force generation models are written in the form

$$\begin{equation} \frac{d\mathbf{z}_{\mathrm{act}}}{dt} = \boldsymbol{\Phi }_{\mathrm{act}}\left(\mathbf{z}_{\mathrm{act}}, [\mathrm{Ca}^{2+}]_{\mathrm{i}}, {SL}, \frac{d\,{SL}}{dt} \right), \qquad T_{\mathrm{a}} = {q}(\mathbf{z}_{\mathrm{act}}), \cssId{texmlid14}{\tag{4.1}} \end{equation}$$

where $\mathbf{z}_{\mathrm{act}}(t) \in \mathbb{R}^{N_{\mathrm{act}}}$ denotes a vector collecting the state variables associated with the dynamics of the contractile proteins and $T_{\mathrm{a}}$ denotes the generated active force. The state variables $\mathbf{z}_{\mathrm{act}}$ may be associated either with the probability of the regulatory proteins of being in a given state Reference 51Reference 70Reference 78Reference 95, or with the elongation of the contractile proteins Reference 72Reference 97, or finally they may be phenomenological variables Reference 51Reference 55. The dynamics of $\mathbf{z}_{\mathrm{act}}(t)$ is driven by that of the calcium concentration ($[\mathrm{Ca}^{2+}]_{\mathrm{i}}$) and of the sarcomere length, denoted by ${SL}$. As a matter of fact, when the sarcomere is elongated (i.e., for large ${SL}$), the generated force tends to increase; moreover, when the sarcomere is shortening (i.e., when ${d\,{SL}}/{dt} < 0$), the generated force is lower than in steady-state conditions (i.e., when ${d\,{SL}}/{dt} = 0$).

In recent years, many efforts have been devoted to the development of force generation models Reference 50Reference 56Reference 70Reference 77Reference 80. The main challenge facing modelers is to describe the complex subcellular processes underlying force generation in a compact way Reference 74Reference 77. In fact, despite several models that describe in a biophysically detailed manner the force generation mechanisms that have been proposed, these are not suitable for numerical simulations because of the overwhelming computational cost associated with their numerical approximation. These difficulties are mainly related to two reasons.

The first reason is the stochastic nature of the processes in play, due to the scale at which they occur (the spatial scales are so small that the effects of thermal agitation cannot be overlooked) Reference 81. For this reason, the output of a force generation model is not a deterministic solution, but rather a probability distribution on the space of the trajectories of the system. Clearly, this dramatically increases the dimensionality of the problem.

The second reason lies in the nearest-neighbor interactions taking place between the proteins that make up the sarcomere Reference 48Reference 77. Because of these interactions, in fact, mean-field approximations, which are widely used in the context of molecular dynamics to lower the computational burden of numerical simulations, cannot be adopted without compromising the validity of the results. In fact, these nearest-neighbor interactions are at the basis of an important phenomenon, namely the cooperativity of contractile units, for which the calcium-activated units favor the activation of neighboring units, thus generating a highly nonlinear relationship between calcium concentration and generated force. These cooperative processes, of fundamental importance for an efficient functioning of the organism, cannot be grasped by mean-field models, in which a single representative unit is considered and the other units are taken into account through the average effect they have on the representative unit Reference 77.

To give a concrete example, let us consider the model proposed in the seminal work Reference 78, in which the authors demonstrated the importance of a spatially explicit description of the sarcomere proteins. In this model, each contractile unit is described by a 4-state Markov chain. Considering a single myofilament with $N = 32$ units, we have $4^N \simeq 2 \cdot 10^{19}$ possible states. Hence, to characterize a probability distribution over this set, we need an order of $10^{19}$ variables. In conclusion, to simulate on a computer the dynamics of this stochastic model, more than $10^5$ petabytes (i.e., $10^{20}$ bytes) would be required just to store the system state in the computer memory. This corresponds to more than 30000 times the storage capacity of the largest supercomputer in the world (up to June 2020).⁠Footnote¹ This situation is not uncommon: a mathematical model is available; however, its practical interest is very limited because of the overwhelming computational cost of its numerical approximation. In these cases, suitable modeling assumptions and/or mathematical tools are employed to derive an approximate, yet computationally feasible, reduced mathematical model.

The TOP500 Project, https://www.top500.org (URL consulted on July 29, 2020).

✖

In the case at hand, in order to capture the cooperative effects without explicitly tracking the joint probability of units of the whole filament, several strategies have been proposed in the literature. Among these, we mention Monte Carlo sampling techniques Reference 96Reference 97 and the reduction of the number of unknowns obtained by grouping together suitable subsets of the states Reference 18Reference 49Reference 95. Alternatively, in Reference 79 the transition rates of the Markov chain were expressed as nonlinear functions of the calcium concentration, to phenomenologically reproduce the cooperative behavior. In Reference 71Reference 72, we introduced a physically motivated assumption of conditional independence of the stochastic processes associated with units that are far from each other along the filament, given the state of the intermediate units. With this assumption, we can describe the protein dynamics through a nonlinear system of nearly 2000 ODEs, thus allowing us to numerically approximate the solution of one heartbeat in just a few seconds of computational time.

A different approach is that of phenomenological models Reference 44Reference 50Reference 55Reference 79: these do not derive equations from first principles, rather they are built by fitting the measured data with simple laws, chosen a priori by the modeler. The numerical solution of these models, typically expressed as systems of few ODEs, feature a lower computational cost than physics-based models. However, they provide a poorer insight than physically detailed models and they are characterized by a hampered predictive power under experimental conditions different from those under which they are built.

5. Modeling cardiac mechanics

During each heartbeat, the heart muscle undergoes large deformations (up to few centimeters). To model the myocardium displacement, the strain of the tissue must be related to the internal stress induced by cardiomyocytes’ contraction and to the pressure exerted by the blood onto the endocardium, the internal surface of the cardiac chambers. The conceptual framework to describe this phenomenon is continuum mechanics, of which we recall basic notions in the following. For more detail, we refer the interested readers to, e.g., Reference 8Reference 58.

Let $\Omega _{0} \subset \mathbb{R}^3$ be an open connected set, which we denote as the reference (or undeformed) configuration, designating the region of space occupied by an elastic body at rest. In our setting this will be the region occupied by the heart muscle, say at the end of the diastole (the time $t=0$ of our simulations). We consider a time-dependent deformation map $\varphi : \Omega _{0} \times [0,{T}] \to \mathbb{R}^d$, such that $\mathbf{x} = \varphi (\mathbf{X}, t)$ (spatial coordinate) represents the position occupied by the point $\mathbf{X} \in \Omega _{0}$ (material coordinate) at time $t$. We thus define the displacement field $\mathbf{d}(\mathbf{X}) \coloneq \varphi (\mathbf{X}) - \mathbf{X}$ and the deformation gradient tensor $\mathbf{F}(\mathbf{X},t) = \nabla \varphi (\mathbf{X},t) = \mathbf{I} + \nabla \mathbf{d}$, where $\mathbf{I}$ denotes the identity tensor and $\nabla$ is the gradient operator in the reference (material) coordinate. The deformation map is assumed to be smooth enough (typically twice continuously differentiable, but weaker regularity is admitted), injective, and orientation preserving (i.e., its Jacobian ${J} = \det \mathbf{F} > 0$ for any $\mathbf{X} \in \Omega _{0}$). By denoting by $\mathrm{Lin}$ the vector space of the linear transformations from $\mathbb{R}^d$ into itself, let us introduce the subset $\mathrm{Lin}^{+}\coloneq \{ \mathbf{A} \in \mathrm{Lin} \text{ s.t.}\ \det \mathbf{A} > 0 \}$.

By Newton’s second law, the displacement field $\mathbf{d}(t)$ satisfies the momentum balance equation,

$$\begin{equation} \rho \, \partial _{tt} \mathbf{d} - \nabla \cdot \mathbf{P} = \mathbf{h} \qquad \text{in } \Omega _{0} \times (0,T], \cssId{texmlid3}{\tag{5.1}} \end{equation}$$

endowed with suitable boundary conditions, where $\rho$ is the mass density of the body and where $\mathbf{h}$ represents an externally applied load (force per unit volume). The tensor $\mathbf{P}$, known as the first Piola–Kirchhoff stress tensor (or simply the Piola stress tensor), encodes the internal stresses of the body. For the sake of model closure, the stress tensor $\mathbf{P}$ is defined via a constitutive law, which is a relationship linking the state of strain of the body with its state of stress. This law can possibly depend on the rate of strain (e.g., in the case of visco-elastic materials). For elastic materials the stress tensor can be written in terms of the strain, i.e., $\mathbf{P} = \mathbf{P}_{\mathbf{F}}(\mathbf{F})$. Here we focus on hyperelastic materials, characterized by a strain energy density $\mathcal{W}$, such that $\int _{\Omega _{0}} \mathcal{W}(\mathbf{X}) dV_{0}$ provides the total elastic energy stored by the body as a consequence of the deformation, where $\mathcal{W}(\mathbf{X}) = \mathcal{W}_{\mathbf{F}}(\mathbf{F}(\mathbf{X}))$ for some $\mathcal{W}_{\mathbf{F}} : \mathrm{Lin}^{+} \to \mathbb{R} \cup \{+\infty \}$. By definition, hyperelastic materials are such that

$$\begin{equation} \mathbf{P} = \frac{\partial \mathcal{W}}{\partial \mathbf{F}}. \cssId{texmlid15}{\tag{5.2}} \end{equation}$$

The passive mechanical response of the heart is significantly anisotropic, due to the presence of fibers. Many anisotropic constitutive laws have been proposed in the literature to describe the cardiac tissue (see, e.g., Reference 39Reference 40Reference 43Reference 92), accounting for the different elastic response along the three directions $\mathbf{f}_{0}$, $\mathbf{s}_{0}$, and $\mathbf{n}_{0}$, the three preferential directions that also characterize the electrical diffusivity of the cardiac tissue (see Section 3). These energies typically feature an exponential dependence on the strain to model the large stiffening of the tissue when it is over-stretched. As an example, the material model of Reference 92 is defined by the strain energy density function

$$\begin{equation} \begin{split} \mathcal{W} &= \frac{C}{2} \left(e^{Q} - 1\right) + \frac{B}{2}({J}-1)\log {J},\\ Q &= b_{ff} E_{ff}^2 + b_{ss} E_{ss}^2 + b_{nn} E_{nn}^2 + 2 b_{fs} E_{fs}^2 + 2 b_{fn} E_{fn}^2 + 2 b_{sn} E_{sn}^2, \end{split} \cssId{texmlid16}{\tag{5.3}} \end{equation}$$

where $C>0$ is a constant and $E_{ab} = \mathbf{E}\,\mathbf{a}_0 \cdot \mathbf{b}_0$, for $a,b \in \{f,s,n\}$, are the entries of the Green–Saint Venant tensor $\mathbf{E} = \tfrac{1}{2}(\mathbf{F}^T\mathbf{F} - \mathbf{I} )$ in the $(\mathbf{f}_{0},\mathbf{s}_{0},\mathbf{n}_{0})$ frame of reference. The constant $B>0$ is the bulk modulus which weights the volumetric term $\tfrac{1}{2}({J}-1)\log {J}$ by penalizing those deformations that would lead to a change of volume occupied by the tissue (${J} = \det \mathbf{F} \neq 1$). The latter term leads to a quasi-incompressible formulation, as little volume variations are allowed. The strongly incompressible formulation is an alternative, in which the balance of the momentum equation Equation 5.1 is coupled with the constraint ${J} = 1$, and the Piola stress tensor is redefined as $\mathbf{P} = \tfrac{\partial \mathcal{W}}{\partial \mathbf{F}} - p \, {J} \, \mathbf{F}^{-T}$, thus yielding a saddle-point problem, wherein the pressure $p$ acts as a Lagrange multiplier.

The cardiac tissue is in fact an active material. This means that its internal stress is not uniquely identified by strain, but rather it can be produced by microscopic mechanisms that turn the chemical energy of ATP into mechanical work. The models presented in Section 4 describe these processes at the microscale and allow to obtain $T_{\mathrm{a}}$, a scalar that represents the magnitude of active force per unit area. This quantity must be related to the macroscopic balance of the momentum equation Equation 5.1, by writing the Piola stress tensor as

$$\begin{equation} \mathbf{P} = \mathbf{P}^{\mathrm{pass}} + \mathbf{P}^{\mathrm{act}}, \cssId{texmlid4}{\tag{5.4}} \end{equation}$$

namely as the sum of a passive term (given by $\mathbf{P}^{\mathrm{pass}} = \tfrac{\partial \mathcal{W}}{\partial \mathbf{F}}$) and an active term $\mathbf{P}^{\mathrm{act}}$ which must be constitutively defined by suitably upscaling the microscopically generated stress. In particular, if we suppose that cardiac muscle fibers, aligned along $\mathbf{f}_{0}$, generate a force per unit area of magnitude $T_{\mathrm{a}}$, directed along the fibers’ direction in the current configuration $\mathbf{f} \coloneq \mathbf{F} \mathbf{f}_{0} / | \mathbf{F} \mathbf{f}_{0}|$, we get

$$\begin{equation} \mathbf{P}^{\mathrm{act}} = T_{\mathrm{a}} \frac{\mathbf{F} \mathbf{f}_{0} \otimes \mathbf{f}_{0}}{| \mathbf{F} \mathbf{f}_{0}|}. \cssId{texmlid17}{\tag{5.5}} \end{equation}$$

Alternative to the active stress approach Equation 5.4 is the active strain one, wherein activation is modeled as a prescribed strain Reference 5Reference 6. We only consider here the active stress approach, which is by far the most widely used in the literature.

For the mathematical well-posedness of nonlinear elasticity models, we refer interested readers to Reference 12Reference 53. Existence results require suitable hypotheses on the energy density function $\mathcal{W}$. The celebrated existence result by J. Ball Reference 12, based on the notion of polyconvexity (a weaker notion than convexity) and on suitable growth and regularity conditions, is based on the direct method of calculus of variations Reference 27. These results, derived for passive materials, cannot be directly applied to heart mechanics, as the cardiac tissue is in fact an active material. The well-posedness of active strain models has been studied in Reference 6, while hypotheses for the well-posedness of active stress models have been derived in Reference 70.

6. Modeling blood circulation

The balance of momentum equation Equation 5.1 must be supplemented with suitable boundary conditions to account for the forces acting between the myocardium and the surrounding environment. In particular, the internal surfaces of cardiac chambers (the endocardium) are in contact with the blood, which exerts pressure. Since the blood dynamics in the chambers is tightly related to that of the whole vascular network, it is necessary to develop models of the entire circulatory system. For this purpose, several blood circulation models, with different degrees of accuracy, have been proposed. These range from three-dimensional fluid-structure interaction (FSI) models, where the blood flow in blood vessels is described by the Navier–Stokes equations Reference 62Reference 65Reference 87Reference 88Reference 94, to zero-dimensional models (whose variables only depend on time, but not on spatial coordinates) Reference 65Reference 75.

In the latter family of models, also known as lumped parameters models, the circulatory system is split into a finite number of compartments and an average pressure is associated with each of them. Then, equations describing the dynamics of these pressure variables and of the blood volume contained in each compartment are derived by the principles of conservation of mass and momentum. Lumped parameters models may be limited to a subset of the circulatory systems Reference 19Reference 36, or may describe the whole closed-loop cardio-circulatory system Reference 16Reference 42Reference 75. In general terms, they are written as systems of ODEs, in the form

$$\begin{equation} \frac{d\mathbf{z}_{\mathrm{circ}}}{dt} = \boldsymbol{\Phi }_{\mathrm{circ}}\left(\mathbf{z}_{\mathrm{circ}}, t\right), \cssId{texmlid5}{\tag{6.1}} \end{equation}$$

where the vector $\mathbf{z}_{\mathrm{circ}}(t) \in \mathbb{R}^{{N_{\mathrm{circ}}}}$ collects the state variables (pressure and volumes). The right-hand side may depend on the time variable, to account for the different phases of the cardiac cycle. A concrete example of Equation 6.1 is available in Reference 75. The zero-dimensional model of blood circulation Equation 6.1 must be suitably coupled with the three-dimensional model of cardiac mechanics (see Section 5). With this aim, we require that the volume enclosed by the ventricular cavity when the domain $\Omega _{0}$ is moved by the displacement $\mathbf{d}$ (which we denote by $\mathcal{V}_{\mathrm{LV}}^{\mathrm{3D}}(\mathbf{d})$) equals the one predicted by the zero-dimensional circulation model (which we denote by $\mathcal{V}_{\mathrm{LV}}^{\mathrm{0D}}(\mathbf{z}_{\mathrm{circ}})$). This coupling is enforced via a Lagrange multiplier, which, in the case at hand, represents the blood pressure inside the ventricular cavity $p_{\mathrm{LV}}$ (more details will be provided in Section 7).

An intermediate class of models for vessel circulation consists of 1D models, in which only the main direction of blood flow is explicitly represented and the equations are averaged in the orthogonal directions. The development of 1D models for blood circulation owes its origin to Leonhard Euler, back in 1775 Reference 30. The mathematical models of the circulatory system also need to account for the dynamics of the four cardiac valves. The latter regulate blood flow between the atria and the corresponding ventricles, between the right ventricle and the pulmonary artery, and between the left ventricle and the ascending aorta. Their sudden opening and closing is driven by the pressure gradient. The motion of cardiac valves can be described with different levels of detail, ranging from three-dimensional FSI models Reference 61Reference 63 to zero-dimensional models Reference 16. For the sake of space, three-dimensional and one-dimensional circulation models will not be further discussed in this paper.

7. A fully coupled cardiac electromechanics model

The mathematical models discussed in Sections 3, 4, 5, and 6 describe different physical processes occurring at different spatio-temporal scales along each heartbeat, which cannot however be understood independently of one another. These are indeed connected though a complex and fascinating web of interactions, feedback loops, and self-regulation mechanisms aimed at preserving the physiological working regime of the heart and responding to external stimuli (see Figure 3). In mathematical terms, in order to describe the whole cardiac function, these models must be coupled to yield a unique system of PDEs and ODEs.

Let us consider for simplicity a single chamber, namely the left ventricle. Hence, the open connected domain $\Omega _{0} \subset \mathbb{R}^3$ represents the region of space occupied by the left ventricle at rest (e.g., at the end of the diastole). Its boundary is split into $\partial \Omega _{0} = \overline{\Gamma }_{0}^{\mathrm{endo}} \cup \overline{\Gamma }_{0}^{\mathrm{epi}}\cup \overline{\Gamma }_{0}^{\mathrm{base}}$, where $\Gamma _{0}^{\mathrm{endo}}$, $\Gamma _{0}^{\mathrm{epi}}$, and $\Gamma _{0}^{\mathrm{base}}$ are disjoint sets respectively denoting the endocardial and epicardial surfaces and the ventricular base, namely the artificial boundary located where the left ventricle geometry is cut (see Figure 3). Then, the fully coupled electromechanical problem consists in finding

$$\begin{equation} \begin{split} & {v} \colon \Omega _{0} \times [0,{T}] \to \mathbb{R} , \quad \mathbf{d} \colon \Omega _{0} \times [0,{T}] \to \mathbb{R}^3 , \quad \mathbf{z}_{\mathrm{ion}} \colon \Omega _{0} \times [0,{T}] \to \mathbb{R}^{N_{\mathrm{ion}}} , \\ & \mathbf{z}_{\mathrm{act}} \colon \Omega _{0} \times [0,{T}] \to \mathbb{R}^{N_{\mathrm{act}}} , \quad \mathbf{z}_{\mathrm{circ}} \colon [0,{T}] \to \mathbb{R}^{{N_{\mathrm{circ}}}} , \quad p_{\mathrm{LV}} \colon [0,{T}] \to \mathbb{R}\\ \end{split} \cssId{texmlid7}{\tag{7.1}} \end{equation}$$

such that we have, in $\Omega _{0} \times (0,T]$,

$$\begin{equation} \left\{ \begin{aligned} &\chi _{m} {C_{m}} \partial _t {v} - \nabla \cdot \left( {J} \mathbf{F}^{-1}\mathbf{D}\mathbf{F}^{-T} \, \nabla {v} \right) + \chi _{m} \mathscr{I}_{\mathrm{ion}}({v}, \mathbf{z}_{\mathrm{ion}}) = \mathcal{I}_{\mathrm{app}}, \\ & \rho \, \partial _{tt} \mathbf{d} - \nabla \cdot \mathbf{P} = \mathbf{0} , \quad \mathbf{P} = \frac{\partial \mathcal{W}}{\partial \mathbf{F}} + T_{\mathrm{a}}(\mathbf{z}_{\mathrm{act}}) \frac{\mathbf{F} \mathbf{f}_{0} \otimes \mathbf{f}_{0}}{| \mathbf{F} \mathbf{f}_{0}|},\\ &\partial _t \mathbf{z}_{\mathrm{ion}} = \boldsymbol{\Phi }_{\mathrm{ion}}({v}, \mathbf{z}_{\mathrm{ion}}), \\ &\partial _t \mathbf{z}_{\mathrm{act}} = \boldsymbol{\Phi }_{\mathrm{act}}\left(\mathbf{z}_{\mathrm{act}}, [\mathrm{Ca}^{2+}]_{\mathrm{i}}(\mathbf{z}_{\mathrm{ion}}), {SL}(\mathbf{d}), \partial _t {SL}(\mathbf{d}) \right) \\ \end{aligned} \right. \cssId{texmlid18}{\tag{7.2}} \end{equation}$$

with boundary conditions

$$\begin{equation} \left\{ \begin{aligned} & \left( {J} \mathbf{F}^{-1}\mathbf{D}\mathbf{F}^{-T} \, \nabla {v} \right) \cdot \mathbf{N} = \mathbf{0} && \, \text{on } \partial \Omega _{0} \times (0,T], \\ & \mathbf{P} \mathbf{N} + K_{\text{epi}} \mathbf{d} + C_{\text{epi}} \partial _t \mathbf{d} = \mathbf{0} && \, \text{on } \Gamma _{0}^{\mathrm{epi}} \times (0,T], \\ & \mathbf{P} \mathbf{N} = -p_{\mathrm{LV}} {J} \mathbf{F}^{-T} \mathbf{N} && \, \text{on } \Gamma _{0}^{\mathrm{endo}} \times (0,T], \\ & \mathbf{P} \mathbf{N} = p_{\mathrm{LV}} \, | {J} \mathbf{F}^{-T} \mathbf{N} | \, \frac{\int _{\Gamma _{0}^{\mathrm{endo}}} J \mathbf{F}^{-T} \mathbf{N} d \Gamma _0}{\int _{\Gamma _{0}^{\mathrm{base}}} | J \mathbf{F}^{-T} \mathbf{N} | d \Gamma _0} && \, \text{on } \Gamma _{0}^{\mathrm{base}} \times (0,T], \\ \end{aligned} \right. \cssId{texmlid11}{\tag{7.3}} \end{equation}$$

and coupled with the following equations in $(0,T]$,

$$\begin{equation} \frac{d \mathbf{z}_{\mathrm{circ}}}{dt} = \boldsymbol{\Phi }_{\mathrm{circ}}\left(\mathbf{z}_{\mathrm{circ}}, p_{\mathrm{LV}}, t\right) , \qquad \mathcal{V}_{\mathrm{LV}}^{\mathrm{3D}}(\mathbf{d}) = \mathcal{V}_{\mathrm{LV}}^{\mathrm{0D}}(\mathbf{z}_{\mathrm{circ}}) \cssId{texmlid19}{\tag{7.4}} \end{equation}$$

with initial conditions (in $\Omega _{0} \times \{0\}$)

$$\begin{equation} \begin{aligned} {v} = {v}_0, \, \ \ \mathbf{d} = \mathbf{d}_0, \, \ \ \partial _t \mathbf{d} = \mathbf{0}, \, \ \ \mathbf{z}_{\mathrm{ion}} = \mathbf{z}_{\mathrm{ion,0}} , \, \ \ \mathbf{z}_{\mathrm{act}} = \mathbf{z}_{\mathrm{act,0}} , \, \ \ \mathbf{z}_{\mathrm{circ}} = \mathbf{z}_{\mathrm{circ,0}}. \\ \end{aligned} \cssId{texmlid8}{\tag{7.5}} \end{equation}$$

We remark that the displacement of the myocardium (encoded in the variable $\mathbf{d}$) affects the conductivity properties of the tissue (thanks to the so-called mechano-electrical feedback), for which Equation 3.3 is rephrased in the deformed configuration, as the tensor

$$\begin{equation} \begin{split} \mathbf{D} &= \sigma _{f} \frac{\mathbf{F}\mathbf{f}_{0} \otimes \mathbf{F}\mathbf{f}_{0}}{|\mathbf{F}\mathbf{f}_{0}|^2} + \sigma _{s} \frac{\mathbf{F}\mathbf{s}_{0} \otimes \mathbf{F}\mathbf{s}_{0}}{|\mathbf{F}\mathbf{s}_{0}|^2} + \sigma _{n} \frac{\mathbf{F}\mathbf{n}_{0} \otimes \mathbf{F}\mathbf{n}_{0}}{|\mathbf{F}\mathbf{n}_{0}|^2}, \end{split} \cssId{texmlid20}{\tag{7.6}} \end{equation}$$

namely the pullback of the conductivity tensor in the reference configuration $\Omega _{0}$ is replaced into Equation 3.2.

The state variable of the ionic model (i.e., $\mathbf{z}_{\mathrm{ion}}$) provides the value of $[\mathrm{Ca}^{2+}]_{\mathrm{i}}$ at each point of the computational domain and at each time $t$, which is an input for the force generation model to be solved at each point in $\Omega _{0}$. To provide the missing input to the latter model, namely the sarcomere length and its time derivative, we observe that the elongation of sarcomeres can be obtained as ${SL}(\mathbf{d}) = {SL}_{0} |\mathbf{F}\mathbf{f}_{0}|$, where ${SL}_{0}$ is the reference sarcomere length.

On its turn, the state variable of the force generation model (i.e., $\mathbf{z}_{\mathrm{act}}$) provides the magnitude of $T_{\mathrm{a}}$, the active force generated at a given point of $\Omega _{0}$ and at a given time $t \in [0,{T}]$, which is employed within the equation describing the mechanics of the myocardium. The associated boundary conditions model the interaction of the cardiac tissue with the surrounding tissue. Specifically, the boundary condition imposed on $\Gamma _{0}^{\mathrm{epi}}$ encodes an elastic ($K_{\text{epi}}$) and viscous ($C_{\text{epi}}$) interaction with the pericardium, a thin membrane enclosing the heart Reference 36. Conversely, the boundary condition enforced on $\Gamma _{0}^{\mathrm{endo}}$ models the stress exerted by the blood contained in the cavity, whose pressure—tracked by the variable $p_{\mathrm{LV}}$—is assumed to be constant in space. The variable $p_{\mathrm{LV}}$ acts in fact as a Lagrange multiplier, enforcing the geometrical compatibility condition $\mathcal{V}_{\mathrm{LV}}^{\mathrm{3D}}(\mathbf{d}) = \mathcal{V}_{\mathrm{LV}}^{\mathrm{0D}}(\mathbf{z}_{\mathrm{circ}})$ at each time $t \in (0,T]$. Finally, the boundary condition imposed on $\Gamma _{0}^{\mathrm{base}}$—known as energy-consistent boundary condition—accounts for the stress exerted by the upper portion of the heart (that we have deliberately disregarded for the sake of simplicity) on $\Omega _0$, consistently with the principles of energy conservation (see Reference 73 for a thorough derivation). We remark that different sets of boundary conditions, representing different means of interaction with the surrounding tissue, could be considered Reference 36Reference 73.

Despite the large amount of work on cardiac electromechanical models, there are still many open questions about the mathematical well-posedness of these systems of equations. In Reference 60 the authors studied the existence of solutions for the fully coupled electromechanical model, considering different kinds of mechano-electrical feedback in an active stress setting. Concerning the active strain setting, instead, in Reference 7, the authors proved existence of weak solutions and uniqueness of regular solutions for a simplified electromechanical model, where the nonlinear elastic model is linearized and where the FitzHugh–Nagumo ionic model Reference 54 is employed. In Reference 15 these results have been generalized to the more realistic Beeler–Reuter Reference 14 and Luo–Rudy Reference 52 ionic models. The latter results are based on the Faedo–Galerkin method and compactness arguments. (A similar approach, named the Galerkin method, stands also at the basis of the numerical approximation and will be outlined in Section 9.)

Part 2. Numerical approximation

8. Why should we turn to numerical approximation

The electromechanical model Equation 7.1–Equation 7.5 is in principle able to encode the EM heart function. As a matter of fact, reality is more complex than that. Indeed, Equation 7.1–Equation 7.5 has been proposed for the left ventricle. A model with a similar mathematical structure can be derived for the right ventricle, as well as for each one of the two atria. The complete EM model would require using different constitutive laws for the atria (their tissue is much thinner than in ventricles), accounting for the electrical signal activation (from the sino-atrial node to the sino-ventricular one, then spreading through the whole myocardium), and coupling the four chambers together; see Reference 65. In this work, we focus only on the left ventricle model Equation 7.1–Equation 7.5, as the associated mathematical considerations will highlight what is relevant for the complete model too.

There are a few other issues that need to be addressed, though. This is precisely where applied mathematicians must put their hands into the “real matter”: How do we generate nonhomogeneous boundary data? How do we provide an initial state to our system (i.e., initial data for our initial boundary value problem)? How do we prescribe a patient-specific problem’s coefficients (like fibers orientation) without subjecting a patient to unnecessary expensive and invasive clinical exams and measurements? Here is where we must depart from the beauty and lightness of PDE theory where problems are invariably formulated under the following incipit: “Let us consider the following PDE with regular enough boundary conditions given by (some Dirichlet or Neumann data) and initial conditions (some suitable initial functions)”. As a matter of fact, boundary and initial functions are often missing for a living system, and even when available they are often incomplete and discrete (that is available only at few selected points) their “regularity” is hard to attribute, actually it might not even be appropriate to talk about functional spaces for boundary or initial data. In practice, this shortcoming (a troublesome issue) is overcome by resorting to suitable mathematical procedures aimed at providing part of the missing data “coherently” with the original physical problem. See for instance Reference 65 and Reference 69. See also Sections 11 and 13 where some of these issues will be addressed in the context of our numerical simulations.

From a theoretical standpoint, the first and foremost question we should address is whether Equation 7.1–Equation 7.5 is well-posed. To our knowledge, a rigorous existence and uniqueness theorem, under realistic assumptions on boundary data, initial conditions, and the problem’s coefficients, is not provable as yet. Partial results do exist concerning the electrophysiology subproblem Equation 3.2 with suitable boundary conditions Reference 34Reference 93, as well as for the passive mechanics problem Equation 5.1 under suitable boundary conditions Reference 12Reference 27. See Reference 7Reference 15Reference 60 for existence results on the coupled electromechanical model under simplifying assumptions. This may raise questions about the actual solvability of problem Equation 7.1–Equation 7.5 and suggests a cautious attitude to mathematicians. On the other side, if we are confident—based on the derivation of Equation 7.1–Equation 7.5 from physics first principles and from biophysically motivated assumptions—that this model represents a faithful description of the reality, then reality suggests to us that “a solution” should exist, and is unique, under physiological (living) conditions. The above heuristic considerations prompt us to try to find numerical solutions (exact, analytic solutions in closed form are impossible to obtain, unless under severe, unrealistic simplistic assumptions).

Generally speaking, the approximation of an initial boundary value problem (IBVP) (either scalar, or vectorial, as the case of Equation 7.1–Equation 7.5) can be operated following different strategies. To avoid unnecessary (and lengthy) generalization, here we will stick with sequential approximations, first in space and next in time (simultaneous space-time approximations can be pursued as well).

9. Space discretization

Space approximation has the aim of reducing the given IBVP to a system of ODEs (i.e., an IVP). This again can be operated according to many different paradigms. Here we will use the Galerkin method that consists in projecting the weak (variational) formulation of our IBVP at every time $t$ into a finite dimensional subspace of the (Sobolev) space, say $V$, to which the exact (weak) solution belongs. More specifically we will use the Galerkin finite element (FE) method for which the finite-dimensional space, say $V_h$, is made of piecewise polynomial functions (either continuous or discontinuous) on a partition of the computational domain into simplexes (triangles in 2D, tetrahedra in 3D); see, e.g., Reference 22Reference 68. The subindex $h$ denotes the characteristic diameter of the simplexes: the smaller the $h$, the finer the FE partition, the higher the dimension, say $N_h$, of the vector space $V_h$, the more accurate (in principle) the numerical solution. Accordingly, we will denote by $y_h(\mathbf{X}, t)$ the FE approximation of the solution $y(\mathbf{X}, t)$ of the IBVP. By suitably defining a basis of the vector space $V_h$, say $\{\varphi _i\}_{i=1}^{N_h}$, the numerical solution is expressed as a linear combination of this basis through $N_h$ (unknown) coefficients that depend on the time variable $t$, that is to say $y_h(\mathbf{X}, t) = \sum _{i=1}^{N_h} \varphi _i(\mathbf{X})\, y_{i,h}(t)$, where