Aorta zero-stress state modeling with T-spline discretization

The image-based arterial geometries used in patient-specific arterial fluid–structure interaction (FSI) computations, such as aorta FSI computations, do not come from the zero-stress state (ZSS) of the artery. We propose a method for estimating the ZSS required in the computations. Our estimate is based on T-spline discretization of the arterial wall and is in the form of integration-point-based ZSS (IPBZSS). The method has two main components. (1) An iterative method, which starts with a calculated initial guess, is used for computing the IPBZSS such that when a given pressure load is applied, the image-based target shape is matched. (2) A method, which is based on the shell model of the artery, is used for calculating the initial guess. The T-spline discretization enables dealing with complex arterial geometries, such as an aorta model with branches, while retaining the desirable features of isogeometric discretization. With higher-order basis functions of the isogeometric discretization, we may be able to achieve a similar level of accuracy as with the linear basis functions, but using larger-size and much fewer elements. In addition, the higher-order basis functions allow representation of more complex shapes within an element. The IPBZSS is a convenient representation of the ZSS because with isogeometric discretization, especially with T-spline discretization, specifying conditions at integration points is more straightforward than imposing conditions on control points. Calculating the initial guess based on the shell model of the artery results in a more realistic initial guess. To show how the new ZSS estimation method performs, we first present 3D test computations with a Y-shaped tube. Then we show a 3D computation where the target geometry is coming from medical image of a human aorta, and we include the branches in our model.

A challenge very specific to patient-specific arterial FSI computations, such as patient-specific aorta FSI computations, is how to use the image-based arterial geometry. The image-based geometry does not come from the zero-stress state (ZSS) of the artery. Special methods targeting cardiovascular MBI and FSI include those designed to account for that. The attempt to find a ZSS for the artery in the FSI computation was first made in a 2007 conference paper [66], where the concept of estimated zero-pressure (EZP) arterial geometry was introduced. The method introduced in [66] for calculating an EZP geometry was also included in a 2008 journal paper on ST arterial FSI methods [32], as "a rudimentary technique" for addressing the issue. It was pointed out in [32,66] that quite often the image-based geometries were used as arterial geometries corresponding to zero blood pressure, and that it would be more realistic to use the imagebased geometry as the arterial geometry corresponding to the time-averaged value of the blood pressure. Given the arterial geometry at the time-averaged pressure value, an estimated arterial geometry corresponding to zero blood pressure needed to be built. Special methods developed to address the issue include the newer EZP versions [20,34,37,38,65] and the prestress technique introduced in [16], which was refined in [18] and presented also in [20,65].
We introduced in [67] a method for estimation of the element-based ZSS (EBZSS) in the context of finite element discretization of the arterial wall. The method has three main components. (1) An iterative method, which starts with a calculated initial guess, is used for computing the EBZSS such that when a given pressure load is applied, the image-based target shape is matched. (2) A method for straight-tube segments is used for computing the EBZSS so that we match the given diameter and longitudinal stretch in the target configuration and the "opening angle." (3) An element-based mapping between the artery and straighttube is extracted from the mapping between the artery and straight-tube segments. This provides the mapping from the arterial configuration to the straight-tube configuration, and from the estimated EBZSS of the straight-tube configuration back to the arterial configuration, to be used as the initial guess for the iterative method that matches the image-based target shape. Test computations with the method were also presented in [67] for straight-tube configurations with single and three layers, and for a curved-tube configuration with single layer. The method was used also in [55] in coronary arterial dynamics computations with medical-image-based time-dependent anatomical models.
In [68], we introduced the version of the EBZSS estimation method with isogeometric wall discretization, using NURBS basis functions. With isogeometric discretization, we can obtain the element-based mapping directly, instead of extracting it from the mapping between the artery and straight-tube segments. That is because all we need for the element-based mapping, including the curvatures, can be obtained within an element. With NURBS basis functions, we may be able to achieve a similar level of accuracy as with the linear basis functions, but using larger-size and much fewer elements, and the NURBS basis functions allow representation of more complex shapes within an element. The 2D test computations with straight-tube configurations presented in [68] showed how the EBZSS estimation method with NURBS discretization works. In [69], which is an expanded, journal version of [68], we also showed how the method can be used in a 3D computation where the target geometry is coming from medical image of a human aorta.
In this article we are introducing a new method for estimating the ZSS. The estimate is based on T-spline discretization of the arterial wall and is in the form of integration-pointbased ZSS (IPBZSS). The method has two main components.
(1) An iterative method, which starts with a calculated initial guess, is used for computing the IPBZSS such that when a given pressure load is applied, the image-based target shape is matched. (2) A method, which is based on the shell model of the artery, is used for calculating the initial guess. The T-spline discretization enables dealing with complex arterial geometries, such as an aorta model with branches, while retaining the desirable features of isogeometric discretization. The IPBZSS is a convenient representation of the ZSS because with isogeometric discretization, especially with T-spline discretization, specifying conditions at integration points is more straightforward than imposing conditions on control points. Calculating the initial guess based on the shell model of the artery results in a more realistic initial guess. To show how the new method for estimating the ZSS performs, we first present 3D test computations with a Y-shaped tube. Then we show a 3D computation where the target geometry is coming from medical image of a human aorta, and we include the branches in our model.
In Sect. 2, we describe the Element-Based Total Lagrangian (EBTL) method, including the EBZSS and IPBZSS concepts. How the initial guess is calculated based on the shell model of the artery is described in Sect. 3. The numerical examples are given in Sect. 4, and the concluding remarks in Sect. 5.

EBTL method
In this section we provide an overview of the EBTL method [67], including the EBZSS concept, and describe the IPBZSS concept and the conversion between the two ZSS.
Let Ω 0 ⊂ R n sd be the material domain of a structure in the ZSS, where n sd is the number of space dimensions, and let Γ 0 be its boundary. Let Ω t ⊂ R n sd , t ∈ (0, T ), be the material domain of the structure in the deformed state, and let Γ t be its boundary. The structural mechanics equations based on the total Lagrangian formulation can be written as Here, y is the displacement, w is the virtual displacement, δE is the variation of the Green-Lagrange strain tensor, S is the second Piola-Kirchhoff stress tensor, ρ 0 is the mass density in the ZSS, f is the body force per unit mass, and h is the external stress vector applied on the subset (Γ t ) h of the boundary Γ t .

EBZSS
In the EBTL method the ZSS is defined with a set of positions X e 0 for each element e. Positions of nodes from different elements mapping to the same node in the mesh do not have to be the same. In the reference state, X REF , all elements are connected by nodes, and we measure the displacement y from that connected state. The implementation of the method is simple. The deformation gradient tensor F is evaluated for each element: where x is the position in the deformed configuration. The deformation gradient tensors for different elements are on different states, but the terms in Eq. (1), including the second term, do not depend on the orientation. Therefore the rest of the process is the same as it is in the total Lagrangian formulation.

IPBZSS
The key idea behind the EBZSS method was that, due to the objectivity, all the quantities seen in Eq. (1) can be computed with any orientation of the ZSS. We can extend the way we see X e 0 to integration-point counterpart of X e 0 . As we did with the EBZSS, we work with the reference domain. With the reference Jacobian Eq. (1) can be rearranged as In our implementation, we use the natural coordinates, with covariant basis vectors where ξ I is the parametric coordinate, and I = 1, . . . , n pd , with n pd being the number of parametric dimensions. The contravariant basis vectors can be calculated with the metric tensor components as where With those vectors, we can express the deformation gradient tensor: and the Cauchy-Green deformation tensor: The Jacobian, J = det F, can be expressed as We can write det C as and from that we obtain We define the covariant basis vectors corresponding to X REF : and the components of the metric tensor are In Eq. (20), we replace g I and g J with and obtain an alternative to the expression given by Eq. (4): The Green-Lagrange strain tensor, where I is the identity tensor, can be expressed with the contravariant basis vectors as The second Piola-Kirchhoff tensor can be expressed with the covariant basis vectors as where S I J can be expressed with the components of the metric tensors. Thus, the inner product δE : S, and all the other quantities, in the weak form given by Eq. (5) can be evaluated without actually using the basis vectors G I . This justifies using (G I J ) k as the integration-point counterpart of X e 0 , with k = 1, . . . , n int , where n int is the number of integration points. Note that G I J is symmetric, and therefore the IPBZSS representation will in 3D have 6×n int parameters for each element.

EBZSS to IPBZSS
Converting the EBZSS representation to IPBZSS representation is straightforward. From given X e 0 we can calculate the covariant basis vectors at each integration point ξ ξ ξ k : and obtain the components of the metric tensor from Eq. (11).

IPBZSS to EBZSS
Converting the IPBZSS representation to EBZSS representation, which we might need for visualization purposes, will, in general, not be exact because the IPBZSS has more parameters than the EBZSS. Given (G I J ) k , we solve a steady-state element-based problem (with f = 0 and h = 0): and the solution to that, in the form X REF + y, will be the EBZSS representation. If the stress calculated from the solution is zero, then the conversion will be exact. We note that to obtain a steady-state solution, we need to preclude translation and rigid-body rotation by imposing 6 appropriate constraints. To do that we first select three control points: A, B and C. We set all three components of y A to be zero and constrain y B to be in the direction The last constraint is y C to be on the plane defined by the vector

An iterative method
Here we assume that we have a reasonably good initial guess for the IPBZSS (see Sect. 3). The iterative method below is used in calculating the IPBZSS that results in the target state associated with the given load. In the iterative method, we estimate F from the ith solution. We simply assume which means The inverse of the deformation gradient tensor from the ith solution can be written as Similarly, the inverse of the target deformation gradient tensor is Thus, we obtain the following equation: Inner-producting both sides of this equation from the right with the covariant basis vectors corresponding to X REF , we obtain which results in The components of the metric tensor are Rearranging the terms, we obtain Thus, we update the components of the metric tensor without actually knowing the basis vectors G I . Once we obtain the IPBZSS at the end of the iterations, we will have the option to convert it to EBZSS and use that in the subsequent y computations.

Initial guess based on the shell model of the artery
An analytical relationship between the ZS and reference states of straight-tube segments was given in [67]. The relationship was called "straight-tube ZSS template" in [69] and was extended to curved tubes. These were for the EBZSS. Here we directly build the IPBZSS instead of building the EBZSS first. This is simpler because with isogeometric discretization, especially with T-spline discretization, specifying conditions at integration points is far more straightforward than imposing conditions on control points. We start with the artery inner surface, which is what the medical images show. Typically, we cannot discern the wall thickness from the medical image. Therefore we first build the inner-surface mesh with T-splines. Then we build a Tspline volume mesh by extruding the surface elements by an estimated thickness.
In our notation here, x will now imply X REF , which is our "target" shape, and X will imply X 0 . We explain the method in the context of one element in the thickness direction. Extending the method to multiple elements is straightforward.

Inner-surface coordinates in the target state
The coordinate system we have here is similar to the one used for the shell modeling in [70]. We note that the "midsurface" of the shell formulation has been shifted to the inner surface here, and • indicates the inner surface. The basis vectors are where α = 1, . . . , n sd − 1, and the third direction is The second fundamental form is defined as and the curvature tensor iŝ Clearly,κ κ κ is symmetric. For a given unit vector t on the surface, we obtain the curvatureκ aŝ If t is a principal direction, andκ is the corresponding principal curvature. The eigenvector can be expressed as Substituting this into Eq. (43) and inner-producting with g α , we obtain where the mixed components indicatê Since the inverse of g γ α exists, κ is an eigenvalue of the matrix defined by the mixed componentsκ α •β , and we will call the two eigenvaluesκ 1 andκ 2 . Forκ 1 >κ 2 we obtain the corresponding eigenvectors: where t 1 and t 2 are unit vectors. From Eq. (43), we writê Sinceκ κ κ is symmetric, Substituting Eqs. (50) and (51) into this, we obtain Thus, the two vectors are orthonormal, and we can express the curvature tensor aŝ κ κ κ =κ 1 t 1 t 1 +κ 2 t 2 t 2 .
Whenκ 1 =κ 2 , an arbitrary orthonormal set of t 1 and t 2 can be used in the above equation. For more details on calculating the eigenvalues and eigenvectors, see Appendix A.

Inner-surface coordinates in the ZSS
Since the principal curvature directions t 1 and t 2 of the target shape are orthogonal to each other, we can build the ZSS shape using those directions. The basis vectors on the inner surface in the ZSS are and the third direction is The stretches corresponding to those directions will beλ 1 andλ 2 . Then the ZSS basis vectors are calculated from That iŝ Because the third direction is orthogonal to t 1 and t 2 , we can reduce these equations tô Substituting Eqs. (48) and (49) into these, we get This can also be written as and from that we calculate the basis vectors as

Wall coordinates in the target state
The position in the target configuration is where 0 ≤ ϑ ≤ h th , and h th is the wall thickness in the target configuration. The basis vectors will vary along the thickness direction as The second and third lines are explained in Appendix 1. The third coordinate is mapped as where −1 ≤ ξ 3 ≤ 1. The basis vector in the third direction is constant as and With that, the components of the metric tensor are

Wall coordinates in the ZSS
The position in the ZSS configuration is where 0 ≤ ϑ 0 ≤ (h th ) 0 , and (h th ) 0 is the wall thickness in the ZSS configuration. The basis vectors will vary along the thickness direction as The curvature tensor in the ZSS configuration iŝ From that, Similar to what we had for t 1 and t 2 , which becomes We substitute Eq. (73) into this and obtain and

Calculating the components of the ZSS metric tensor at each integration point
For an integration point ξ ξ ξ , we can obtain the components of the metric tensor as The third coordinate can be obtained from Assuming incompressible material, J = 1, where The components of the matrix tensors are given by Eqs. (74) and (89).

Design of the ZSS
The design parameters are the principal curvatures κ 0 1 and κ 0 2 , and the stretchesλ 1 andλ 2 for each principal curvature direction. Those parameters can be determined fromκ 1 and κ 2 of the target configuration.
As proposed in [69], the two principal directions are seen as circumferential and longitudinal directions, andκ 1 is in the circumferential direction, giving us Here φ is the opening angle, which is seen after a longitudinal cut, based on artery experimental data [71]. The stretch in that direction,λ 1 , corresponds to λ θ in [69] and that is deter- Thickness (mm) mined from the 2D computations in [69]. We assume that in the longitudinal direction the ZSS configuration has zero curvature, κ 0 2 = 0. The stretch in that direction,λ 2 , corresponds to λ z in [69]. If at an integration pointκ 2 ≈κ 1 , we set κ 0 2 = κ 0 1 andλ 2 =λ 1 . That makes the assignment of the principal directions less consequential.
where We assume the pressure associated with the target shape is 92 mm Hg. The initial guess for the iterations is determined as described in Sect. 3, with φ = 5 2 π andλ 2 = 1.05. After the iterations, explained in Sect. 2.5, for comparison purposes, we convert the IPBZSS representation to EBZSS representation, with the method described in Sect. 2.4. With the EBZSS,

Y-shaped tube
The target state of the Y-shaped tube is shown in Fig. 1. The end diameters of the tube are 20, 14 and 10 mm. Figure 2 shows the T-spline mesh. The mesh is based on a mixture of cubic and quartic T-splines. The wall thickness distribution is smooth, outcome of solving the Laplace's equation over the inner surface, with Dirichlet boundary conditions at the tube ends, where the value specified is 0.1 times the end diameter. Figure 3 shows the thickness distribution. The volume mesh is built with one element (cubic Bézier element) in the thickness direction. The number of control points and elements are 5,180 and 2,592. Since the IPBZSS cannot be visualized, we show it using the EBZSS representation. Figure 4 shows the initial guess for the IPBZSS. Figure 5 shows an element in the target state and the corresponding IPBZSS initial guess.
We iterate and obtain the converged IPBZSS. Figure 6 shows y computed from that. The maximum value of y is 1.279×10 −14 mm. Figure 7 shows the von Mises strain, computed from the IPBZSS initial guess and from the converged IPBZSS. The converged IPBZSS element is shown in Fig. 8. Figure 8 shows that the opening angle and the longitudinal stretch do not change much between the IPBZSS initial guess and the converted IPBZSS. However, the circumferential stretches are somewhat different. The initial

Remark 1
Thickness (mm) Fig. 13 Patient-specific aorta geometry. Wall thickness distribution guess for the circumferential stretch was based on the 2D computations reported in [69]. Alternatively, we can estimate the stretch by an analytical solution, similar to the one described in [70].
We also compute y after converting the converged IPBZSS to EBZSS. Figure 9 shows y computed that way. The maximum value of y is 1.626×10 −2 mm. Figure 10 shows the von Mises strain. There is no visible difference between the strains obtained from the IPBZSS directly and after conversion to EBZSS.

Patient-specific aorta geometry
The target state of the patient-specific geometry is shown in Fig. 11. Figure 12 shows the T-spline mesh.
The wall thickness distribution is smooth, outcome of solving the Laplace's equation over the inner surface, with Dirichlet boundary conditions at the tube ends, where the value specified is 0.08 times the end diameter. In some parts of the branched area the thickness exceeds the radius of curvature, and there we reduce the thickness to 0.8 times the radius of curvature. Figure 13 shows the thickness distribution.
The volume mesh is built again with one element (cubic Bézier element) in the thickness direction. The number of control points and elements are 9,244 and 4,360. Figure 14 shows the initial guess for the IPBZSS. We again iterate and obtain the converged IPBZSS. Figure 15  shows y computed from that. The maximum value of y is 1.163×10 −13 mm. Figure 16 shows the von Mises strain, computed from the IPBZSS initial guess and from the converged IPBZSS. We again compute y after converting the converged IPBZSS to EBZSS. Figure 17 shows y computed that way. The maximum value of y is 1.057×10 −1 mm. Figure 18 shows the von Mises strain. These results indicate that the EBZSS obtained by conversion from the IPBZSS is also a reasonable representation of the aorta ZSS.

Concluding remarks
We have introduced a new method for estimating the ZSS required in patient-specific arterial FSI computations, where the image-based arterial geometries do not come from the ZSS of the artery. The estimate is based on T-spline discretization of the arterial wall and is in the form of IPBZSS. The method has two main components. (1) An iterative method, which starts with a calculated initial guess, is used for computing the IPBZSS such that when a given pressure load is applied, the image-based target shape is matched.
(2) A method, which is based on the shell model of the artery, is used for calculating the initial guess. The T-spline discretization enables dealing with complex arterial geometries, such as an aorta model with branches, while retaining the desirable features of isogeometric discretization. The desirable features of higher-order basis functions of isogeometric discretization include being able to achieve a similar level of accuracy as with the linear basis functions, but using larger-size and much fewer elements, and being able to represent more complex shapes within an element. The IPBZSS is a convenient representation of the ZSS because with isogeometric discretization, especially with T-spline discretization, specifying conditions at integration points is more straightforward than imposing conditions on control points. Calculating the initial guess based on the shell model of the artery results in a more realistic initial guess. To show how the new ZSS estimation method performs, we first presented 3D test computations with a Y-shaped tube. Then we presented a 3D computation where the target geometry was coming from medical image of a human aorta, and the model included the aorta branches.