A Mathematical Spline-Based Model of Cardiac Left Ventricle Anatomy and Morphology †

: Computer simulation of normal and diseased human heart activity requires a 3D anatomical model of the myocardium, including myoﬁbers. For clinical applications, such a model has to be constructed based on routine methods of cardiac visualization, such as sonography. Symmetrical models are shown to be too rigid, so an analytical non-symmetrical model with enough ﬂexibility is necessary. Based on previously-made anatomical models of the left ventricle, we propose a new, much more ﬂexible spline-based analytical model. The model is fully described and veriﬁed against DT-MRI data. We show a way to construct it on the basis of sonography data. To use this model in further physiological simulations, we propose a numerical method to utilize ﬁnite differences in solving the reaction-diffusion problem together with an example of scroll wave dynamics simulation.


Introduction
Integrative mathematical models of complicated hierarchic physiological systems, such as the heart, allow the use of a computational mathematical approach to describe these systems from the molecular level to the macro level, which characterizes the structure and functions of the whole heart and/or certain heart chambers, including the LV.In the framework of numerical experiments, a comprehensive analysis of such ventricular models provides a way to elucidate the mechanisms of normal and pathological heart performance; these models also provide a way to predict potential methods of effective therapy, such as correcting mechanical and electrical cardiac function disturbances.
Preferential paths for the propagation of an electrical excitation wave in the human ventricular myocardium are associated with muscle fibers in tissue.The speed of an excitation wave along a fiber is several times higher than that across the direction of a fiber.Thus, 3D simulation of electrophysiological and mechanical activity of the heart requires an anatomical model that includes a field of myofiber directions.
Recently, several models of the electrical and/or mechanical activity of the whole heart or its chambers have been proposed [1][2][3][4][5][6][7][8][9][10][11][12].The most valuable of these are based on a detailed description of cardiac anatomy and the fiber orientation field, which is crucial for a correct representation of the physiological function of the heart.

Theoretical Models of Cardiac Anatomy
The theoretical approach to the modeling of heart chamber, e.g., the LV, architectonics arises from the attempt to clarify a general principle or governing rule that determines the orientation of the fibers in the ventricle wall.Such theoretical models may be especially useful when studying basic mechanisms of normal and pathologic heart performance.
The discovery of rotational anisotropy in the myocardium was one of the first fundamental findings in qualitative heart anatomy.The tangents to the myofibers turn, if we look at them transmurally, and the turning angle equals approximately 120 • .Streeter compared this picture with a Japanese fan [20].Torrent-Guasp suggested that we consider the myocardium of both ventricles as a "flattened rope" or a band that has two endings and makes a few turns [28].According to Streeter [20], myofibers are geodesic lines on a set of nested toroids, and outer toroids are more prolate, while inner ones are more round.Besides this model, there are many models concerning a few myocardial layers (usually three, namely superficial, middle and deep; see, for example, [29]).One can find a more detailed description of the history of cardiac qualitative anatomy in [30].
One group of models constructs the architectonics of myocardial fibers by using the analogy of muscular skeletal fibers.Different models include from one to four systems of myocardial bundles [31][32][33][34][35].However, some researchers point to restrictions of the approach [36][37][38].
One of the key questions concerning the LV myocardium is related to its layers, namely their existence, and if they exist, their number and positions.As a rule, models with sheets have a number of layers inserted one into another [29].Nevertheless, a layers-based model exists in which the number of sheets is indefinite and the LV is not divided into internal and external parts.The model in question is from an early work [28], published in 1973, in which Torrent-Guasp proposed a hypothesis on how cardiac fiber orientation could be reproduced using geometrical transformations.However, this idea has never been formulated in mathematical terms.
Let us now discuss mathematical models of the heart.Theoretical approaches to cardiac architectonics include rule-based methods [44][45][46] and the mapping of fibers from animal hearts to human heart geometry [47].One of the most recent rule-based methods is a Laplace-Dirichlet algorithm [46], which takes a noisy DTI-derived fiber orientation field as input data and yields, firstly, the transmural and apicobasal directions for the entire myocardium and, secondly, a smooth and continuous fiber orientation field.Another approach was used by Peskin, who derived the fiber orientation field from the principles of mechanical equilibrium [48].
Chadwick considered a cylindrical LV model and specified the helix angle linearly depending on the position of the point in the LV wall [49].A spheroidal LV model, which involved the linear dependence of the helix angle on the distance between the point and the endocardium, was proposed by Beyar et al. in [2].
Another theoretical approach was created by Arts et al. in 1992 [50].Their model was an ellipsoidal LV with a piecewise quadratic law for helix angle change.The mechanical adaptational principle suggested by Arts et al. in 1982 [51] was used to quantify the orientation of muscle fibers via the helix fiber angle distribution.

Relation to Other Anatomical Ventricular Models
Torrent-Guasp in [28] suggested that the LV can be represented as a set of embedded surfaces, i.e., muscular layers, covered by curves representing the myofibers.This idea was illustrated using Pettigrew's proposal [52] of these surfaces having the form of a planar semicircle with a number of curves drawn on the semicircle.The sheets thereafter were transformed into conical surfaces.Finally, the heart was considered as a family of such conical surfaces inserted into each other.It may be noted that Streeter mentioned the importance of these ideas, but did not formalize them.
The present work develops two analytical LV models based on Pettigrew's anatomical concept.One of the models is axisymmetric [53].The other model is built on power functions and is not symmetric [54].Heart anatomy was described in these works using a family of inserted non-conical surfaces, which allowed us to represent the smooth geometry of the LV properly.A family of lines on these surfaces was used to represent cardiac fibers.The entire model was expressed in terms of analytical functions.
Our experience shows that in many situations, a well-fitting non-symmetric LV model cannot be made using the approach from [54].We faced this when we tried to make "power-function" models using sonography and tomography data of human hearts.The limitation of the approach is that the LV model cannot fit shapes of some real hearts accurately enough.
The new model uses splines instead of power functions and, as we will show in this paper, is flexible enough to fit real medical visualization data.At the same time, the new model uses analytical calculations, so some of the derivatives necessary for solving reaction-diffusion systems can be found exactly.In the proposed model, both the anatomy and fiber direction field are defined analytically.We calculate fiber slope angles in a local coordinate system and compare them with in vitro experimental data on human hearts.
In the next section, we describe the details of our model construction.In Sections 3 and 4, we compare its output with experimental measurements of anisotropy in the human and canine LV.In Section 5, we propose a technique for constructing a model based on sonography data, which is especially important for clinical applications.Sections 6 and 7 provide a numerical method for simulating the electrophysiological activity of the LV and an example of the usage of the numerical method.Finally, Section 8 is a discussion on the results obtained.

Construction of the LV Model
In order to define the LV form, we use a special coordinate system (γ, ψ, φ), where the variable γ ∈ [0, 1] corresponds to the position of a point in the LV wall layer: γ = 0 is the endocardium; γ = 1 is the epicardium; ψ ∈ [0, π/2] is an analogue of latitude; ψ = 0 is the upper plane part of the LV model (fibrous ring and valve zone); ψ = π/2 is the LV apex; and φ ∈ [0, 2π) is an analogue of longitude.
Initial data for the model construction are: 1. meridians φ i , i = 0, 1, . . .n − 1, where measurements are taken; 2. coordinates ρ epi i,j , z epi i,j , i = 0, 1, . . .n − 1, j = 0, 1, . . .n epi i − 1, of marked points on the epicardium; 3. coordinates ρ endo i,j , z endo i,j , i = 0, 1, . . .n − 1, j = 0, 1, . . .n endo i − 1, of marked points on the endocardium.Now, we describe our algorithm to construct the LV model.First, for each meridian of the epicardium, we connect the marked points by an interpolation curve (Figure 1).We obtain a set of n curves in this step.Second, for each latitude ψ, we connect the n points on the curves by a closed interpolation curve.Algebraically, this means that we use a periodic spline.We obtain an epicardial surface in this step.Next, we do the same for the endocardium and, finally, we fill the space between the epi-and endocardium by a linear function of γ.
This schema can be formalized using the following equations.To transform the special coordinates into cylindrical ones (ρ, φ, z), we use a simple formula for z-coordinate transformation: where Z is the LV height and h is the LV wall thickness at the apex.Using this formula, we calculate ψ-coordinates of the marked points: As mentioned above, γ = 0 on the epicardium, and γ = 1 on the endocardium.
In each section φ i on the epicardium (endocardium), we make an epicardial (endocardial) curve: Here, spline is a cubic spline with a zero second derivative at the endpoints.The notation spline(x; x i → y i ) means a spline constructed by two given arrays of data, x i as an argument and y i as a response variable.Then, we interpolate to obtain a surface between the curves: Here, we use a periodic cubic spline.Finally, we fill the entire LV wall: Usually, if n meridional sections are considered, the angles between them are uniform, so φ i = iπ/n, i = 0, 1, . . .n − 1.

Spiral Surfaces
Hierarchically, the LV model is subdivided into muscle layers.In its turn, each layer is filled by myofibers.All of the layers have a spiral-like shape, so we call them "spiral surfaces" (SS).These surfaces have the same equation in special coordinates: where φ max is the SS twist angle (which is the same for all SS), and different values of φ min ∈ [0, 2π) determine different surfaces.The equation of the SS in the cylindrical coordinates is shown below (see Equations ( 1) and ( 5)).

Filling a Spiral Surface with Fibers
We have utilized Pettigrew's proposal [52] as we did in our previous models [53,54].This proposal has been chosen because it can be formalized easily, and it describes both myofibers and sheets.Myocardial fibers are obtained as images of chords Y = const, Y ∈ [0, 1) of sector P ≤ 1, Φ ∈ [πγ 0 , πγ 1 ] (the chords are parallel to the diameter Y = 0) on the SS (Figure 2).Here, γ 0 and γ 1 (0 ≤ γ 0 < γ 1 ≤ 1) are parameters necessary to fit the sub-epicardial and sub-endocardial fiber angles.The parameter of any chord is the polar angle Φ ∈ [Φ 0 , Φ 1 ], where Φ 0 = max(arcsin Y, πγ 0 ) and Φ 1 = min(π − arcsin Y, πγ 1 ).The mapping of a chord point (P, Φ) to an SS point is defined by the formulae below (Figure 3):  For instance, if γ 0 = 0 and γ 1 = 1, the semicircle diameter is transformed into a fiber that begins on the basal epicardium, descends to the apex (Φ = π/2), then ascends and ends on the basal endocardium.Images of shorter chords are located closer to the LV base and have shorter lengths.
An algorithm for calculating the fibre direction vector in points of the model is given in Appendix A.

Fitting the LV Form
In [54], a fitting of the more rigid model to a real canine heart dataset was done with satisfactory results.Yet, fitting the dataset of a human heart into that model was less successful, partially due to the limitations of the model's parameters.In this study, the same dataset of the human heart was used for the flexible model.The dataset is accessible online at [55].
The fitting procedure began with finding the LV axis Oz.Then, we sectioned the LV by N = 12 meridional half-planes φ i = 2πi/N, i = 0, 1 . . .N − 1, passing this axis, and we manually marked points on the epi-and endocardium in each section.The mean numbers of points were 10 for the sub-endocardium and 11 for the sub-epicardium.In total, taking into account the coincidence of two apical points in all of the sections, there were 227 points.After setting the points, we used periodic cubic φ-splines and cubic ψ-splines.The results obtained are shown in Figures 4 and 5.

Methods for Model and Experimental Data Comparison
In the work [54], the theoretical model and experimental data were compared along normals to epicardial meridional sections as described in [20].However, this method is not useful: normals do not have to intersect the endocardium, and they may intersect each other, so the task of finding a normal passing through a point may result in multiple solutions.In the present article, we use the special coordinates to construct straight pinning lines that do not intersect each other and always intersect the endocardium.These lines have a simple equation: ψ = const, φ = const.
We will now describe the comparison procedure.To compare angles along such a pinning line, one needs to specify a point A on the epicardium.Let its special coordinates be γ = 0, ψ = ψ * and φ = φ * .Let us consider a corresponding meridional section φ = φ * of the model, semi-plane Π.The line lies in Π and intersects the endocardium at a point B, having coordinates γ = 1, ψ = ψ * and φ = φ * .On the segment AB, we set k equidistant points, including its ends, so that A = A 1 , A 2 , . . . ,A k = B.The position of a point A i on the segment AB is defined by the variable: (for the endocardium x = 0; for the epicardium x = 1).We have to draw an SS through every point The problem of finding such an SS is reduced to solving the SS Equation ( 6) with respect to φ min : Strictly speaking, there can be no points from the tomogram exactly on semi-plane Π; therefore, we selected points lying no further than ∆ = 1 mm from the straight line AB and inside the dihedral angle |φ − φ * | ≤ ∆ φ = 0.1 rad= 5.7 • .Figure 6.The definition of the local coordinate system is shown on the (left) of the image.Oxyz is the global Cartesian coordinate system.The blue axis is normal, u, to the epicardium.The red axis is parallel tangent, v.The dark green axis is a meridian tangent, w.The colorful surface is the epicardium; the color depends on altitude z.The curve inside is a model fiber, and it does not lie on the epicardium.The normal axis intersects the curve at a point.On the (right), the definitions of the true fiber angle α and the helix angle α 1 are shown.The thick, dashed line is a tangent to a myofiber segment constructed at the origin of the coordinates.The dashed-and-dotted lines are projections of the myofiber tangent.(Copied from [53] (Figure 9) with a deletion of angle α 2 .) Streeter proposed specifying fiber direction using a local coordinate system (u, v, w) and two angles, "true fiber angle" α and "helix angle" α 1 [20] (p.78) (see Figure 6); these angles are sufficient for specifying the direction of a fiber at a particular point.The true fiber angle reflects fiber direction totally, with both helix and transverse components, while the helix angle neglects the transverse part.The axis u is a normal to the epicardium pointed away from the LV; w is a meridian, i.e., an epicardial tangent lying in a meridional semi-plane and pointed upwards, and v is a parallel, i.e., vector w × u.Angle α ∈ [0, π/2] is between a fiber and the parallel, and angle α 1 ∈ [−π/2, π/2] is between the fiber projection on the plane uv and the parallel.
We compared the two angular characteristics of fiber directions with the experimental data for a human heart in two meridians (one meridian lies in the LV free wall; another lies in the IVS) in the upper, middle and lower parts of the LV wall.

Results of Comparing the Model with DT-MRI Data
In this section, we compare the directions and angles of myocardial fibers in the model and DT-MRI data of one human and one canine heart.

Comparison of the Model with Human Heart Data along Straight Pinning Lines
The following parameter values, which are common for all meridians, were used: LV height Z = 68 mm, LV wall thickness at the apex h = 8 mm, SS twisting angle φ max = 3π, endocardial fiber parameter γ 0 = 0.13 and epicardial fiber parameter γ 1 = 0.9.
Graphs of the dependency of angles α, α 1 on a point position on pinning lines are shown in Figures 7-12.We will now analyze the results obtained.At the upper and middle LV areas (e.g., Figures 7A and 8A), one can see that the vertical axis does not go through the center of the horizontal LV sections, but it is situated closer to the IVS.The axis is positioned there because it must go through the apical LV area, and the LV apex projection to its basal plane is not situated at the center of the base.If one moves the axis to the base center, then the apex is far from the axis in one of the meridional sections; therefore, we cannot fit the LV wall shape using the central position of the axis.
Let us consider the fiber slope angles in one of the LV free-wall meridians.In the upper LV part (see Figure 7), the true fiber angle α (Panel C) in the model accurately reproduces the DT-MRI data.Angle α descends from 75 • on the endocardium to approximately 15 • at the middle of the wall, then it reaches 65 • on the epicardium.The helix angle (Panel D) in the model is also quite close to the experimental data.The middle part of the free wall, which is determined by height (see Figure 8), shows an essentially large dispersion of the both angles' values in the subendocardial zone (x < 0.2).These angles' values are similar to those in the basal zone, so the model can reproduce them well.The angles at the lower part of the LV free wall (see Figure 9) are predicted by the model with slightly lesser accuracy.
In the upper part of the IVS (see Figure 10), the model behaves in a different way than the experimental data.In the middle part of the IVS (see Figure 11), the model reproduces the angles reasonably accurately.Figure 12 shows that the model simulates the angles in the interior part of the myocardium better than in the exterior part.
The transmural variation of angle α is parabolic both in the DT-MRI data and in the model.This angle is locally maximal at the subendo-and subepicardium, and it is minimal, but not zero at the middle of the wall, where the helix angle is zero.Angle α is not zero there due to the traverse fiber angle.
The dashed black lines in Figures 7-12 demonstrate the results of another modeling approach described in [46].We see that the functions are linear with respect to x, and the angles do not depend on the region.Bayer et al.'s model has zero transmural angle (traverse angle α 3 in Streeter's notation), so its angle α = |α 1 |.
If one considers the fibers in the radial direction, our model (like the model from [53]) imitates the distinctive arrangement of fibers in the ventricular wall (see Figure 14 of [53]).As we mentioned above, this arrangement was called the "Japanese fan" by Streeter (see Figure 42C of [20]).

3D Verification
We compared the fiber directions of the human tomography dataset and the model in all DT-MRI points located inside the model LV wall, and we made a histogram of the results (Figure 13).The mean angle between the two directions is 21.1 • , and the median is 15.2 • .The standard deviation is 18 • .Furthermore, we used our model to fit a canine DT-MRI dataset, which is accessible through the link [55] given in Section 2.3.The following parameter values, which are common for all meridians, were used: LV height Z = 90 mm, LV wall thickness at the apex h = 12 mm, SS twisting angle φ max = 3π, endocardial fiber parameter γ 0 = 0.18 and epicardial fiber parameter γ 1 = 0.9.A 3D comparison of the dataset and the model was performed.This comparison showed the following results: mean 19.5 • , median 16.1 • and standard deviation 14.1 • , which is very close to the human heart comparison.
The degree of error in the angles consists of two parts.One of them is the modeling error caused by inaccuracies in the model.This error can be seen in Figures 7-12, where the blue curves are out of the red point clouds.The second part is caused by errors of measurement and by the dispersion of data.Even if a blue curve is in the middle of a point cloud, such errors accrue.To identify these parts of errors, we partly averaged the true fiber angle (α) and helix angle (α 1 ) data using the Smoothn function in MATLAB [56], so we obtained data without outliers.The averaging is controlled by the smoothing parameter s ≥ 0. If s = 0, the algorithm makes no changes in the input data.The noise in the data diminishes with the growth of s.Full averaging was impossible because an increase in the smoothing parameter above s ≈ 5 led to an oversmoothing without an essential noise reduction.We used s = 5, and the results of the 3D comparison are shown in Tables 1 and 2. There, "st.dev." means "standard deviation".We can conclude that the partial averaging of the input data decreases the error sufficiently.The mean error in the fiber angles is about 10-15 • for the data without outliers.

Constructing a Model Based on Sonography Data
The proposed model of the LV form can be constructed based on sonography data.We made one such measurement using data from two-and four-chamber views of one patient's heart.In Figure 14, we show a comparison of two models: the old model based on power functions and the present spline-based one.The model was made using 44 marked points, from which 22 were endocardial and 22 were epicardial.Since the posterior (right on the image) part of the LV wall has a complex curved shape, it cannot be approximated accurately by the first model, but it can be well approximated by the second one.The anterior (left on the image) part of the LV wall has a simpler shape, so both models can be fitted into it well.

Numerical Method for Solving Reaction-Diffusion Systems on the Model
Electrophysiological processes in the myocardium are usually simulated using reaction-diffusion systems of partial differential equations.In the monodomain case, they are: Here, u is the transmembrane potential; v is the vector of other phase variables (its components differ in different models); the diffusion matrix D has elements D ij = D 2 δ ij + (D 1 − D 2 )w i ( r)w j ( r), i, j = 1, 2, 3; and w is the unit vector of the fiber direction.
The boundary condition of zero flux of the potential through the myocardial boundary looks like: where n is a normal vector to the boundary.In the case of complex cardial shape, implementation of the boundary condition is difficult, and it requires either a special method of Laplacian computation near the boundary in Cartesian coordinates or a special coordinate system where the boundaries are coordinate surfaces.
The LV base, epicardium and endocardium are just coordinate surfaces in the special coordinates linked with this model.Indeed, the base has the equation ψ = 0, while the equation γ = const describes the epicardium and endocardium.For the symmetrical LV model, we proposed an algorithm [58], which has been used in a study of scroll wave dynamics [59].For the previous (power-functions based) non-symmetrical model [54], a modified algorithm has been described elsewhere [60].That algorithm can be used in the present case as well without any changes.However, since it was not published in English, we depict its dissimilarity with the symmetrical case method briefly here.Formulae for the Laplacian and boundary conditions in the special coordinates are given in Appendices B and C.

The Method to Rarefy the Computation Mesh
Let us consider a mesh node with indices i, j, k of special coordinates γ, ψ, φ.In the case of the symmetric LV model, we propose that its γ-and ψ-neighbors have indices i ± 1, j ± 1.Its φ-neighbors have indices k ± Q[j], where Q is an array of natural numbers, and the differences are found so that the Cartesian distances between nodes are within some given limits.In the case of the non-symmetric LV model, we have to introduce six arrays instead of one, Q dir var , where var ∈ {γ, ψ, φ}, dir ∈ {+1, −1}.The meaning of Q dir var is that the index of the γ-neighbor of the node in direction The values of Q have to be found once, at the initialization stage.If for a node (i, j, k), at least one of Q dir var [i][j][k] cannot be found, this node is marked "out-of-mesh".We compute the potential in this node not by a finite difference method, but by interpolation, using the node's closest "in-mesh" neighbors.
The flowchart in Figure 15 depicts all of the steps of construction and the use of the model in electrophysiological simulations.
Mark points on the epiand endocardium.
Convert their cylindrical coordinates into special ones.Connect points on each meridional curve by ψ-splines.Connect all the epicardial curves using periodic φ-splines.
Connect all the endocardial curves using periodic φ-splines.
Calculate fibre directions.Find curves and surfaces, if necessary.Correct the parameters, if necessary.
Make mesh(es) for the numerical method and visualisation.
Set electrophysiological and other parameters.Calculate electrophysiological activity.
For each section, do: Set an axis and a flat base.Find h, Z.
Set N meridional sections.

An Example of the Spline-Based Model's Practical Usage in Electrophysiological Simulations
The model was used to study the dynamics of scroll waves of electrical excitation.Such waves are solutions of classical reaction-diffusion systems, and they correspond to dangerous arrhythmias, such as paroxysmal ventricular tachycardia and ventricular fibrillation.Our simulation was based on the Aliev-Panfilov model [61] of electrical excitation in the myocardium and the aforementioned human heart DT-MRI dataset and anatomical model (see Sections 2. 3 and 4).The electrophysiological model is dimensionless and has the following form: where u is the transmembrane potential, v is the conductivity for the K current in the membrane, D 1,2 are diffusion coefficients along and across fibers and k, a, 0 , µ 1 , µ 2 are cell model parameters.These parameters had the values D 1 = 7.12 mm 2 /ms, D 2 = D 1 /9, k = 8, a = 0.1, 0 = 0.01, µ 1 = 0.12, µ 2 = 0.3.A boundary condition provided zero flux through the LV boundary.Initial conditions were ; and u = 0, v = 0 (resting state) elsewhere.We used the spatial mesh size dr = 0.5 mm; the time step dt = 0.01 ms; and the dimension coefficients 5.4 ms, as the model time unit, and 6.2 mm, as the model length unit, based on the period and wavelength of a spiral wave simulated on a flat isotropic surface.The duration of the simulation was 10 s.
The results we obtained are shown in Figures 16 and 17.From these results, we can see that the solitary filament drifted towards the base (ψ = 0), and then, several filaments appeared.The number of filaments slowly increased from 1 to 9 and then decreased to 5-6.This observation means that a tachycardia paroxysm transformed into a fibrillation.

Discussion
This section analyses the algorithm, compares it with other models and techniques, examines its methods of verification and concerns the usage and further development of the constructed model.

Advantages and Disadvantages of the Proposed Model
The model adequately reproduces the fiber angles in the LV free wall and in the middle zone of the IVS of the human heart.Nevertheless, the data agreement in the apical zone of the IVS is only qualitative; in the IVS upper part of the human heart, the model yields results that differ from the experimental data.
Our model has the following merits: • It has relatively few parameters for shape, fibers and sheets.

•
The LV apex is smooth thanks to the special choice of function z.This is a useful feature for integration methods sensitive to the smoothness of the boundary.

•
It uses only simple one-variable splines, and no two-or three-parametric splines are required.

•
The fitting of shape and fiber directions are independent tasks in this model.

•
It is flexible enough to fit not only normal, but pathological LVs.

•
It yields not only fibers, but also sheets.
Limitations of the model include the following: • Similar to our previous models, its base is flat, which is not the case in real mammal anatomy.

•
The fibers are not geodesic lines in the model sheets.It is difficult to determine if this is a disadvantage of the fibers, the sheets or both.

•
Model fibers end at the base.However, this demerit can be amended by combining this model with the toroids-based one described in [62].

•
It is unclear whether the model is generalizable for both ventricles.

The Present Model Compared with Other Qualitative Models
The ventricular myocardium's structure and its qualitative description have been explored for several centuries.Reviews of experimental findings and theoretical conceptions can be found, for example in [23,63] and other reviews in the supplementary volume [64].The present model is conceptually close to Krehl's nested toroids to the extent that the LV myocardium, except the papillary muscles, is filled by non-overlapping surfaces, and each surface is filled by non-overlapping curves.The distinctive feature of the proposed model is the spiral shape of the surfaces.Whether the SSs are close enough to real heart sheets is a topic for future research.
In [20], Streeter proposed the geodesic principle of winding fibers on the sheets.This principle holds in the toroidal model, but not in the present one.However, if we keep the fibers only, the surfaces may be adjusted to them.This interesting problem remains for future research: How should we draw a set of non-overlapping surfaces on the given set of curves so that the curves are geodesic on the surfaces?Some qualitative models of the myocardium describe its structure on the levels of cells, myofibers and myolaminae (see, for example, [22,29,65]); such models examine the interconnections and arrangement of cardiomyocytes and other types of cardiac cells.However, the main task of our model is to yield fiber and sheet normal directions in all points of the LV myocardium, but not to represent the amazingly complex micro-and macro-structure of the cardiac muscle.A union of our model with a structural tissue-level model, which includes Y-junctions of the cells, capillaries, fibrous tissue and other elements, would be of great practical value.
A key feature of the proposed model is its numerical, quantitative nature, so it should be compared with other quantitative models and data.

Comparison with Experimental Data and Quantitative Models
Some experimental techniques that can be used to verify our model have been mentioned in our previous work [54].These techniques include DT-MRI (already used here), as well as micro-CT [66], histological investigation [67] and quantitative polarized light microscopy [30].These methods provide fiber orientation data that may be used for the construction of anatomical computational models in different ways.

•
Rule-based models, i.e., the models formed on the basis of some constitutive rule [44][45][46][47], can be compared against the data, which allows them to be verified.
A comparison with a recent rule-based method, the Laplace-Dirichlet algorithm [46], and with a wrapping-based myocardium model [70] was done in our previous paper [54].The idea of using splines to fit cardiac imaging data is not new.An example of such LV shape model, constructed with B-splines, together with a review of spline-based models, can be found in [71].A common weak point of many such models, however, is a lack of computing fiber direction field.
A quantitative comparison with a model from [46] can be made by measuring the mean difference angle between the experimental and model data.Bayer et al. in [46] report that their result is 23 • ± 20 • .Our result on the human heart dataset is 21 • ± 18 • , and on the canine dataset, it is 19.5 • ± 14.1 • , which is comparable with their achievement.
Our model fits the shape of the LV and has only three parameters (γ 0 , γ 1 , φ max ) to adapt fiber directions.Fiber directions can also be fitted through widely-used fitting functions with many more parameters, as was done in [72].In that work, fiber direction field in mice was approximated in spheroidal coordinates by a product of two algebraic and one trigonometric polynomials.A certain advantage of that approach is its flexibility.The reported root-mean-square error was about 2.3 • -10 • .Yet, that method usually requires about 10 parameters for each algebraic and 10-40 for trigonometric polynomials, whereas our technique uses only three parameters.Furthermore, the parameters in [72] have unclear geometrical or physical meaning, unlike ours.Moreover, the numbers of those parameters for different hearts vary, so the coefficients cannot be averaged to obtain recommended species-level parameter sets.We hope that our parameters can be averaged or fitted on the species level or at least for the norm and pathologies of particular species.This will allow researchers to construct fiber direction fields based on the shape of patients' hearts without DT-MRI data.
Let us discuss the level of model error in the fiber directions.How can local and global cardiac characteristics be affected if the fiber angles change?
In an article by Campbell et al. [73], the electromechanics of dogs' hearts was modeled using different cell-and tissue-level parameter values.One of the authors' simulations included three models with different values of the total transmural fiber angle α 1 , i.e., 112 • , 145 • and 170 • , while all other parameters were the same.The mean transmural difference of the helix angle between the middle model and the extreme models was about 12 • .The authors discovered that such a fiber angle distribution affects the range of fiber strains and weakly changes hemodynamics and transmural delay in activation.
In a conference work by Wong et al. [74], the sensitivity of the average displacement of some points of the heart to epi-and endocardial fiber angles was studied.The angles were changed at increments of 20 • from 20 • -80 • .The analysis showed that a minimal change of the angle (though 20 • is an essential angle) leads to a significant variation in the displacement of the points.Our model has special parameters for fitting the epi-and endocardial fiber angles, so such a significant average error in these angles is unlikely to appear.In Sánchez's PhD thesis [75], linear and cubic laws of the helix fiber angle dependence on wall position were used in electromechanical simulations.The maximal difference between the angles was 23 • , and it led to significant changes in maximal contraction time and total activation time.
Baron et al.'s conference paper [76] is focused on the effect of fiber orientation on LV stroke volume.A change of 10 • in epi-and endocardial helix angles leads to the following effects: an 18% difference in the atrioventricular plane displacement, quite a small difference (≈ 2 mm Hg) in blood pressure in the LV, a 5% difference in the minimal LV-cavity volume, a 3.5% difference in the ejection fraction and a 2% difference in pressure-volume work.
Summarizing all of the aforementioned works, we conclude that fiber (helix) angle differences or errors on the order of 10 • can sufficiently change the local parameters of cardiac activity, but affect its global characteristics lightly.
The goal of our model is to predict the average field of fiber directions for a given species.Therefore, the model error has to lie in the intraspecies limits.For example, the limits reported in [73] for dogs (12 • ) are close to the error achieved by our model (10.9 • ; see Table 1), which justifies the relevance of our technique.

Further Development and Usage of the Model
The model does not include the uppermost part of the ventricular myocardium, the basal ring, where myofibers form torus-like shapes (Figure 14 of [28]).Moreover, some of the data available on the apical region show that myofibers there make toroidal layers, similar to basal ones (Figure 6 of [28]).
This model has been constructed based on and verified against only two DT-MRI datasets, one canine and one human.Such validation is limited, so its development is an important topic of our future research.
The spline apparatus enables researchers to use this model in cases of complex wall geometry, such as in patients with old myocardial infarctions, and other situations where the wall has loci of thickening and thinning.
The analytical description of cardiac geometry has been used in developing a new numerical method for the study of the electrophysiological activity of the LV.The model can also be utilized to generate different anisotropic properties of the heart, to alter the LV shape (by changing the model parameters) and to study their influence on cardiac electrical and mechanical functions.

Appendix B. Formulae for the Laplacian in the Special Coordinates
Let us denote the special coordinates (γ, ψ, φ) uniformly as (ξ 0 , ξ 1 , ξ 2 ) and consider matrices: where v = (v 0 , v 1 , v 2 ) is the unit fiber direction vector.These matrices are linked by the relations: The Laplacian can be written as: where: and q kl are elements of matrix Q: where n is the normal vector to the LV surface, can be written in the form: n T DJ γ ∂u ∂γ + n T DJ ψ ∂u ∂ψ + n T DJ ϕ ∂u ∂ϕ = 0, (C2) where J γ,ψ,ϕ are the vector-columns of derivatives of the special coordinates with respect to Cartesian ones: and so on.We use the method of fictitious nodes and make an additional layer of nodes outside the LV model.Values of the potential are found based on first-order derivatives with respect to the special coordinates.This method allows us to calculate the Laplacian in all non-fictitious nodes, be it points on the LV boundary or in its depth, on a consistent basis.
For the endo-and epicardium, Equation (C2) has to be solved for ∂u ∂γ .Then, we find the potential at the fictitious node behind the endocardium (in the LV cavity) or the epicardium, respectively.
For the LV base, we transform Equation (C1) to form: The required derivatives can be found numerically.The fictitious node is located above the LV base.

Figure 1 .
Figure 1.Vertical (meridional) sections of the LV free wall (on the left) and the IVS (on the right) of a human heart.The points represent DT-MRI data; the solid red line represents the model epicardium; the dashed blue line represents the model endocardium.The circles show the interpolation (marked) points given by a user.The horizontal axis is ρ; the vertical axis is z.On the right panel, we see a papillar muscle in the RV cavity (vertical one, ρ = 25, . . ., 45 mm) and an RV free wall (inclined).

Figure 4 .
Figure 4. Epicardial and endocardial surfaces of the model that was constructed based on DT-MRI data.(left) The epicardium is opaque; (right) The epicardium is semi-transparent.Color shows z-coordinates.

Figure 5 .
Figure 5.A spiral surface (SS) with fibers on it.Two views are shown.The left of the figure depicts a view from the top and side of the SS, while the right of the image provides a side view with the mid-myocardial part on the left and the epicardial layer on the right.Color shows z-coordinates.

Figure 7 .
Figure 7. Fiber angles in the model and in the experimental data.The basal area (ψ = 10 • ) of the LV free wall is shown.(A) A horizontal LV section.The points are myocardial points from a human heart DT-MRI scan.(B) A meridional LV section.The solid (dashed) curve is the model epicardium (endocardium), and the points are myocardial points from the DT-MRI scan.(C,D) The angles α and α 1 , respectively.The X-axis displays the position of points in the wall depth; 0 corresponds to the endocardium, and 1 corresponds to the epicardium.The points show the experimental data, while the solid curves show our model data, and the dashed curves display results yielded by the model from[46].

Figure 8 .
Figure 8. Fiber angles in the model and experimental data.The LV free wall in the middle area (ψ = 25 • ) is shown.The conventional signs are the same as in Figure 7.

Figure 9 .
Figure 9. Fiber angles in the model and experimental data.The LV free wall in the apical area (ψ = 60 • ) is shown.The conventional signs are the same as in Figure 7.

Figure 10 .
Figure 10.Fiber angles in the model and experimental data.The IVS in the basal area (ψ = 5 • ) is shown.The conventional signs are the same as in Figure 7.

Figure 11 .
Figure 11.Fiber angles in the model and experimental data.The IVS in the middle area (ψ = 30 • ) is shown.The conventional signs are the same as in Figure 7.

Figure 12 .
Figure 12.Fiber angles in the model and experimental data.The IVS in the apical area (ψ = 65 • ) is shown.The conventional signs are the same as in Figure 7.

Figure 13 .
Figure 13.Histogram of angles between fibers in experiment and model data.The human heart dataset is used.X-axis is the angle; Y-axis is the number of points in the segments; the segment [0, 90 • ] was divided into 100 segments of equal length.

Figure 14 .
Figure 14.Ultrasound imaging as a base for model construction.(left) A sonographic two-chamber view of a patient's heart and two measurements: the first for the power-function-based model (white lines) and the second for the spline-based model (blue lines show the endocardium; red lines show the epicardium).The patient suffered from an old myocardial infarction, and a scar is present.The marked points are shown as crosses.This view gives data for two model meridians.The model axis and base are shown as grey lines.The image has been rotated so that the model axis Oz is vertical.(right) A spline-based model of this patient's LV.The view is semitransparent to show both the epicardial and endocardial surfaces.The surface mesh was made using Ani3D-AFT software [57].

Figure 15 .
Figure 15.Flowchart displaying the steps for the construction and use of the spline-based LV model.

Figure 16 .
Figure 16.The special coordinates ψ (top) and φ (bottom) of the scroll wave filaments.

Table 1 .
Statistics of the error in the true fiber angle.

Table 2 .
Statistics of the error in the helix angle.