Application of Hyperelastic Nodal Force Method to Evaluation of Aortic Valve Cusps Coaptation: Thin Shell vs. Membrane Formulations

: Coaptation characteristics are crucial in an assessment of the competence of reconstructed aortic valves. Shell or membrane formulations can be used to model the valve cusps coaptation. In this paper we compare both formulations in terms of their coaptation characteristics for the ﬁrst time. Our numerical thin shell model is based on a combination of the hyperelastic nodal forces method and the rotation-free ﬁnite elements. The shell model is veriﬁed on several popular benchmarks for thin-shell analysis. The relative error with respect to reference solutions does not exceed 1–2%. We apply our numerical shell and membrane formulations to model the closure of an idealized aortic valve varying hyperelasticity models and their shear moduli. The coaptation characteristics become almost insensitive to elastic potentials and sensitive to bending stiffness, which reduces the coaptation zone.


Introduction
The human aortic valve prevents regurgitation (backward blood flow), i.e., coapts during diastole. The functionality of the aortic valve decreases with age due to valve sclerosis (calcification), and heart valve diseases are called the "next cardiac epidemic" [1].
One of the most appealing approaches to treating aortic valve sclerosis is through a reconstruction of the aortic valve using chemically treated auto-pericardium [2]. This procedure is low cost, has no immune response and provides minimal degradation of new cusps (leaflets), according to post-operation clinical follow-up.
The key problem for any aortic valve reconstruction procedure is the optimal size and design of new cusps: reliable valve coaptation in its diastolic state contradicts the undesirability of oversized leaflets. The state-of-art design is mostly based on geometrical models [3] or expert decisions [4] that do not account for the anatomical patient-specific characteristics and mechanical properties of the cusp material. An optimal design based on personalized mathematical models is in demand for planing surgery at the preoperative stage and, therefore, has restrictions regarding computational time.
There is a vast amount of the literature devoted to the mathematical modeling of native aortic valve functionality [5] using various models, such as fluid-structure interaction (FSI) models, structural models and simplified models. The most accurate are FSI models; however, they are prohibitively expensive for surgery-planning due to their computational complexity. For instance, an FSI simulation of the aortic valve dynamics during one cardiac cycle may take 195 h [6].
The competence of a reconstructed aortic valve is provided by its coaptation characteristics in the diastolic state. Therefore, the cusps design optimization should be based on coaptation characteristics, which do not depend on the utilized models, FSI or dry structural models ignoring the blood flow [6].
Our approach is to estimate coaptation characteristics (area, profile, heights) in the diastolic state by finding the quasi-static equilibrium of the aortic valve closed by the diastolic aortic pressure. A similar approach is suggested for intraoperative use by surgeons to test postrepair valve competence under physiological pressures [7].
Structural dry-models are still computationally expensive for surgery planning. For instance, the simulation of static systolic mitral valve closure takes 73-98 min on 16 processors [8]. Simplified mass-spring models (MSM) have been proposed for mitral valve surgical planning [8,9], and aortic valve surgical planning [10]. MSM computation is very fast and the coaptation zones obtained by FE (finite elements) and MSM are very similar for the aortic valve [10]. However, for the mitral valve, MSM moderately overestimates the displacement of the leaflets towards the atrium as compared to ground truth data [8]. The main disadvantage of the MSM approach is that it is hard to ensure that you are dealing with an elastic material of interest. In particular, the choice of spring stiffnesses for MSM is unclear in the case of hyperelastic materials. Nevertheless, MSM is believed to be a good trade-off between simulation accuracy and time-efficiency [8].
In our previous work [11], we proposed the use of a hyperelastic nodal force (HNF) method instead of MSM. This method belongs to the class of P 1 finite element discretizations and the easy implementation of various hyperelastic materials. Since the HNF method was applied to the membrane formulation, its computation time was small enough. Moreover, the HNF method was verified on deformations of nonlinear membranes [12]. To account for the bending stiffness of chemically treated pericardium [13], we will consider a pericardium-made leaflet as a thin shell. The compromise between accuracy and computational complexity is based on combination of rotation-free finite elements describing bending stiffness [14] and the HNF method describing membrane behaviour. The rotationfree finite elements were proposed to approximate the curvature of a triangulated surface on a patch defined by a central triangle and its three adjacent triangles [15]. In the present paper, we study how bending stiffness affects the coaptation characteristics of an idealized aortic valve made of human pericardium.
The outline of the paper is as follows. In Section 2, we introduce main notations and describe our method. In Section 3, we verify the method on popular benchmarks and study the coaptation characteristics of an idealized aortic valve. In Section 4, we address the limitations of our study and discuss future work.

Shell Kinematics
We consider the deformation of a thin hyperelastic shell and employ the Kirchhoff-Love assumption that straight lines normal to the reference mid-surface remain normal to the mid-surface after deformation [16].
Let mapping Ψ : A → Ω 0 ∈ R 3 define the mid-surface for an initial (before deformation) configuration. Here A ⊂ R 2 is a compact set. A shell at the initial configuration S 0 is defined as where H is an initial thickness of the shell and N is the unit normal to the mid-surface at the initial configuration. Analogously, we can define an actual (deformed) configuration S t of the shell where ψ : A → Ω t ∈ R 3 defines the mid-surface at the deformed configuration, λ = h/H relates the thickness at the undeformed and deformed configurations, h is the thickness of the deformed shell, and n is the unit normal to the mid-surface at the deformed configuration.
Convective coordinate systems for the initial (undeformed) and the current (deformed) configurations are: respectively. Hereafter, we denote the partial derivative ∂ f /∂ξ α by f ,α . The geometry of the mid-surface is given by the metric tensor with components and the curvature tensor with components The deformation gradient F for the shell is given by where In the theory of thin shells, the right Cauchy-Green deformation tensor C = F T F = ∑ 3 i,j=1 g i · g j G i ⊗ G j is reduced to the strain measure [15,16] which means Thus, tensor C is split into the surface (in-plane) part C S = ∑ 2 α,β=1 c αβ G α ⊗ G β and the out-of-plane part C N = λ 2 N ⊗ N.

Constitutive Equations
The shell is assumed to be made of a hyperelastic material. This means there is an elastic potential (strain energy density function) W(C) such that the second Piola-Kirchhoff stress tensor S has the form [17]: The elastic potential of an isotropic material is a function of three invariants [17]: For incompressible hyperelastic shells I (3D) 3 = 1, and we rewrite the 3D invariants in terms of the surface invariants I 1 , J [18]: therefore, According to (9), the constitutive equation of an incompressible isotropic hyperelastic material is written in the form of the plane stress state where (A) αβ denotes αβ-entry of matrix A.

Weak Formulation
We consider equilibrium of a thin shell subject to mixed boundary conditions and external forces with density b. Let the boundary of the mid-surface ∂Ω t be split into two parts, where u = x − X is the displacement of the mid-surface, n t is the unit outward normal to ∂Ω t , T = (1/ det F) F S F T is the Cauchy stress tensor,ū andt are given displacement and tension on corresponding boundaries. The equilibrium equations are based on the virtual work principle which states that the virtual work of the external forces δW ext is equal to the change in the total stored strain energy δU [19]: Taking (8) and (13) into account, we can rewrite (15)-(17) as follows The first term in (18) represents in-plane membrane behavior, whereas the second term characterizes the bending part (bending stiffness of the thin shell). The nodal hyperelastic force method [12] is applied to the discretization of the first term, and the bending part is discretized by the approach proposed by Oñate et al. [15].

Discretization
Given a consistent triangular mesh of the initial mid-surface configuration Ω 0 , we apply a combination of the P 1 finite elements for the membrane part and the rotation free bending elements for the shell bending part to the approximate solution of Equations (15), (17) and (18).
Let deformation of a triangle T P with vertices P 1 , P 2 , P 3 into a triangle T Q with vertices Q 1 , Q 2 , Q 3 be defined via mapping x(X) (see Figure 1a). Denote the area of the undeformed triangle T P by A P , the area of the deformed triangle T Q by A Q , and the unit normal to the plane of triangle T Q by n. The Jacobian of the deformation J = A Q /A P is one of the surface invariants.

Discretization of the Membrane Part
The hyperelastic nodal forces method is the P 1 finite element discretization of (15)- (17) in the membrane formulation [12].
The membrane part of the deformation of T P is discretized by where u h is a P 1 finite element displacement and the corresponding nodal force for the i-th node of triangle T P is Here, the elastic potential W h | ξ 3 =0 is constant over the triangle T P since it is computed from the linear displacement u h . Formula (23) stems from the one-point quadrature for finding l ij in (19).

Bending Part
The bending part of the deformation of T P is discretized by The contribution of the bending part to the total nodal force is obtained by rotationfree triangular shell elements [15]. The curvature tensor is recovered from nodal displacements of a patch Π P composed of the central element (T P ) and three adjacent elements (in Figure 1b triangles 1, 2, 3). In notations of Figure 1b the local axis ξ 1 is directed along vector t 1 connecting nodes 1 and 2, the local axis ξ 3 is directed along the unit normal n to the plane of T P , the local axis ξ 2 is directed along n × t 1 .
The constant curvature field on each triangle T P is computed by where A P is the area of the central element T P of the patch Π P in the undeformed configuration, ψ h ,αβ is the second partial derivative of P 2 vector function ψ h with given nodal values Q i in the patch Π P (Enhanced Basic Shell Triangle, EBST) [15].
The numerical integration in (25) gives the following expression for the curvature tensor components where denotes the outer normal to the k-th edge of the central element T P in the initial configuration (see Figure 1b), l k and E k denote the length and the middle point of the k-th edge, respectively. Equation (26) relates the curvature variation δκ with vector δQ = (n · δQ 1 , n · δQ 2 , n · δQ 3 , n · δQ 4 , n · δQ 5 , n · δQ 6 ) T and the curvature matrix B b [15] Bending forces at the i-th node of the patch Π P are calculated by which implies the patch Π P contribution to the nodal bending forces

Discretized Equilibrium Equations
The equilibrium Equations (15)-(21) form a local nonlinear system for new positions Q i of mesh nodes. Assembling local contributions of all triangles T P results in the static equilibrium for the i-th node where Σ i is the set of triangles sharing the i-th node and Π i is the set of patches containing the i-th node. Equation (30) form the system of algebraic equations . Parameter τ < 1 is chosen to avoid the divergence of iterations, in all experiments τ = 0.5 except the conservative loading of a cantilever where τ = 0.9.

Numerical Results
We begin the study of our formulation with several benchmark problems for thin shells [20].

Cantilever Subjected to End Shear Force
We consider a cantilever made of the St-Venant-Kirchhoff (SVK) material subjected to the end shear force 0 ≤ P ≤ P max . Geometry setting and elastic constants (Young modulus E and Poisson's coefficient ν) are presented in Figure 2. For this benchmark ν = 0, and we cannot use (11) to express the elastic potential in terms of the surface invariants (I 1 , J), since (11) was derived under the incompressibility condition ν = 1/2. For arbitrary ν, the plane stress state s 13 = s 23 = s 33 = 0 implies the following elastic potential for SVK material in terms of (I 1 , J): The clamped boundary condition is implemented following [15]. If a triangle T P is adjacent to a free edge of Ω 0 , its patch Π P is incomplete and we estimate the curvature tensor from the extra condition of vanishing linearized moments Here, n f = (n 1 , n 2 ) T and t f = (−n 2 , n 1 ) T denote normal and tangent vectors to the free edge in the initial configuration Ω 0 , M (0) denotes a 2 × 2 matrix with components The Equation (33) stems from (20), (8) and linearization of the second Piola-Kirchoff stress tensor S.
We consider conservative and non-conservative loadings of the cantilever [21]. The first one is a standard load P with the fixed direction normal to the initial mid-surface Ω 0 . The second one is a load P perpendicular to the current mid-surface Ω t .
The equilibrium is searched on three quasi-uniform unstructured triangular meshes with mesh sizes h 1 = 0.1, h 2 = h 1 /2, h 3 = h 1 /4. The load-deflection curves for the conservative and non-conservative loadings are shown in Figures 3 and 4, respectively. The curves converge to the reference curves [20,21] as the mesh is refined, with a relative pointwise error of less than 2%.

Slit Annular Plate Subjected to Lifting Line Force
The next benchmark addresses a slit annular plate made of SVK material [20]. The line force P is applied at one end of the slit, while the other end of the slit is fully clamped, see Figure 5. The free edge and clamped boundary conditions were imposed similarly to the cantilever benchmark. The equilibrium is searched on three quasi-uniform unstructured triangular meshes with mesh sizes h 1 = 0.2, h 2 = h 1 /2, h 3 = h 1 /4. The load-deflection curves are shown in Figure 6. The relative pointwise error with respect to the reference curve [20] is less than 2% as well.

Aortic Valve in Diastolic State
The competence of a reconstructed aortic valve is given by cusps coaptation characteristics. The cusps (leaflets) can not interpenetrate, they interact with each other forming coaptation zones. The coaptation zone and other coaptation characteristics for a cusp are schematically represented in Figure 7.  To estimate how bending stiffness affects a reconstructed tricuspid aortic valve competence, we consider idealized geometries of three leaflets and the aortic root [10]. Each leaflet is represented by a semicircle with 20 mm diameter. The leaflet is fixed along its semi-circumference to the inside wall of a cylinder (aortic root) with 60 mm circumference. The idealized aortic valve model combines three identical semicircular leaflets, arranged circumferentially around the inside wall of the cylinder. Such a valve model defines the initial configuration. Each leaflet is subjected to diastolic blood pressure of 80 mm Hg. The boundary conditions are presented in Figure 8. The numerical solution is computed on two quasi-uniform unstructured triangular meshes with mesh sizes h 1 = 0.5 mm (coarse), We compare the coaptation characteristics of the idealized aortic valve obtained by the membrane and the thin shell models. The cusps are assumed to be cut from treated human pericardium. The mechanical properties of fresh and treated human pericardium are not clear: according to [22], fresh and treated human pericardium is isotropic whereas animal pericardium is anisotropic; according to [23], fresh and treated human pericardium is anisotropic. Glutaraldehyde-treated autopericardium is known to be significantly stiffer than native autopericardium. In the present work, we compare only isotropic incompressible materials (St.Venant-Kirchhoff, neo-Hookean, Gent), and vary their elastic moduli. For the last two materials, the elastic potentials are where µ = E/3 is the shear modulus, J m is a material constant. Models (31), (34), (35) benefit from the physical meanings of their parameters. Each valve cusp is represented by an oriented triangulated surface, whose position in the 3D space provides the static equilibrium of the three closed cusps under blood diastolic pressure. The equilibrium implies the following for each free mesh node with index i where F p i , F s i , F c i are the force of diastolic blood pressure, the elastic force and the contact force due to reaction of the other leaflets [11], respectively.
To find the static equilibrium, we apply Algorithm 1 to our inhouse code. Contact detection and contact forces are computed in lines 11-31 of the Algorithm, following the idea of real-time simulation [24] and Soft-Soft contact processing algorithm [25].

Algorithm 1 Solution of the static equilibrium equation and computation of the contact forces
Require: initial positions of leaflets meshes nodes x 0 1: set relaxation coefficient λ, contact distance parameter d c 2: set maximal allowed number of relaxation iterations maxits, convergence parameter ε abs 3: set i = 0 4: repeat 5: compute elastic forces F e ← F e (x i ) 6: compute pressure forces 8: if m > d and v · n > 0 then 26: 27: modify F c at node N:  Table 1 demonstrate that bending stiffness reduces significantly the coaptation zone. Comparison of the coaptation characteristics computed on the coarse and the fine grids indicates that the mesh convergence is achieved.
Stiff cusps (E = 10 MPa) do not coapt when bending stiffness is incorporated. Though they do not coapt properly even in the membrane formulation due to improper leaflet design. Less stiff cusps provide the coaptation. The coaptation height depends on elastic modulus E but does not depend on the material model which confirms the same observation in [11]. In Table 2 we present the CPU time for different models and materials on the coarse mesh. The simulations were run on a laptop Intel Core with i5-8250U CPU 1.60 GHz. For the analysis of the computational work we refer to the next section.  Figure 9a,b demonstrates the closed aortic valve for the neo-Hookean material with E = 1 MPa; the other materials give similar shapes. Convexity at the boundary and the incomplete coaptation for the shell-based model near the commissure points (Figure 9c) are caused by the clamped boundary condition. One can also observe the incomplete coaptation at the center of the leaflets for the shell-based model, due to te insufficient length of the cusps. The membrane-based model for the materials with E = 1 MPa provides the complete coaptation; the gap between the cusps is due to visualization of the mid-surfaces, which are distanced from each other at 0.5 mm, the cusps thickness.

Discussion
Coaptation characteristics are crucial in assessment of competence of reconstructed aortic valves. The valve cusps coaptation can be evaluated via shell [26] or membrane [10] formulations. In this paper, we compare both formulations in terms of coaptation characteristics for the first time. In particular, we study the impact of bending stiffness on the aortic valve closure. Our numerical thin shell model is based on a combination of the hyperelastic nodal forces method and the rotation free finite elements. The numerical shell model is verified on several popular benchmarks for thin shell analysis. The relative error with respect to the reference solutions does not exceed 1-2%.
We apply our numerical shell and membrane formulations to model the closure of the idealized aortic valve varying hyperelastic models and their shear moduli. The coaptation characteristics become sensitive to the bending stiffness and almost insensitive to the elastic potential.
Our study has several limitations. First, the idealized leaflet design [10] leads to incomplete coaptation for both membrane and shell models of stiff cusps. Advanced leaflet designs [3,4], result in better coaptation characteristics. Second, our isotropic constitutive models may overestimate the impact of bending stiffness. Third, the clamped boundary condition may require modifications in order to avoid regurgitation in the vicinity of the commissures. Fourth, the CPU time of computing coaptation for the shell formulation is prohibitively large for patient-specific leaflet shape optimization and our numerical method for shells has to be more efficient in terms of CPU time.
In our future work, we shall address different known leaflet designs, anisotropic materials, more realistic boundary conditions accounting leaflet suturing and employ patient-specific geometry of the aortic root based on segmentation algorithms [27].