Medical-image-based aorta modeling with zero-stress-state estimation

Because the medical-image-based geometries used in patient-specific arterial fluid–structure interaction computations do not come from the zero-stress state (ZSS) of the artery, we need to estimate the ZSS required in the computations. The task becomes even more challenging for arteries with complex geometries, such as the aorta. In a method we introduced earlier the estimate is based on T-spline discretization of the arterial wall and is in the form of integration-point-based ZSS (IPBZSS). 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 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. The method has two main components. 1. An iteration technique, which starts with a calculated ZSS initial guess, is used for computing the IPBZSS such that when a given pressure load is applied, the medical-image-based target shape is matched. 2. A design procedure, which is based on the Kirchhoff–Love shell model of the artery, is used for calculating the ZSS initial guess. Here we increase the scope and robustness of the method by introducing a new design procedure for the ZSS initial guess. The new design procedure has two features. (a) An IPB shell-like coordinate system, which increases the scope of the design to general parametrization in the computational space. (b) Analytical solution of the force equilibrium in the normal direction, based on the Kirchhoff–Love shell model, which places proper constraints on the design parameters. This increases the estimation accuracy, which in turn increases the robustness of the iterations and the convergence speed. To show how the new design procedure for the ZSS initial guess performs, we first present 3D test computations with a straight tube and a Y-shaped tube. Then we present a 3D computation where the target geometry is coming from medical image of a human aorta, and we include the branches in the model.

A challenge specific to patient-specific arterial FSI computations is how to use the medical-image-based arterial geometry, which does not come from the zero-stress state (ZSS) of the artery. The special methods targeting cardiovascular MBI and FSI include those intended to account for that. The task becomes even more challenging for arteries with complex geometries, such as the aorta. 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 [31], as "a rudimentary technique" for addressing the issue. It was pointed out in [31,66] that quite often the medical-image-based geometries were used as arterial geometries corresponding to zero blood pressure, and that it would be more realistic to use the medical-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. The special methods developed to address the issue include the newer EZP versions [20,33,36,37,64] and the prestress technique introduced in [16], which was refined in [18] and presented also in [20,64].
A method for estimation of the element-based ZSS (EBZSS) was introduced in [67] in the context of finite element discretization of the arterial wall. The method has three main components. 1. An iteration technique, which starts with a calculated ZSS initial guess, is used for computing the EBZSS such that when a given pressure load is applied, the medical-image-based target shape is matched. 2. A technique 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 straight-tube 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 ZSS initial guess for the iteration technique that matches the medical-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 [54] in coronary arterial dynamics computations with medical-image-based time-dependent anatomical models.
The version of the EBZSS estimation method with isogeometric wall discretization, using NURBS basis functions, was introduced in [68]. 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. Because all we need for the elementbased 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 fewer elements, and the NURBS basis functions enable 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], how the method can be used in a 3D computation where the target geometry is coming from medical image of a human aorta was also shown.
In the method introduced in [70], the estimate is based on T-spline discretization of the arterial wall and is in the form of integration-point-based ZSS (IPBZSS). 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. The method has two main components. 1. An iteration technique, which starts with a calculated ZSS initial guess, is used for computing the IPBZSS such that when a given pressure load is applied, the medical-image-based target shape is matched. 2. A design procedure, which is based on the Kirchhoff-Love shell model of the artery, is used for calculating the ZSS initial guess.
In this article we increase the scope and robustness of the method by introducing a new design procedure for the ZSS initial guess. The new design procedure has two features. (a) An IPB shell-like coordinate system, which increases the scope of the design to general parametrization in the computational space. (b) Analytical solution of the force equilibrium in the normal direction, based on the Kirchhoff-Love shell model, which places proper constraints on the design parameters. This increases the estimation accuracy, which in turn increases the robustness of the iterations and the convergence speed. To show how the new design procedure for the ZSS initial guess performs, we first present 3D test computations with a straight tube and a Y-shaped tube. Then we present a 3D computation where the target geometry is coming from medical image of a human aorta, and we include the branches in the model.
In Sect. 2, from [70], we describe the Element-Based Total Lagrangian (EBTL) method, including the EBZSS and IPBZSS concepts. The IPB shell-like coordinate system and the related mesh generation are described in Sect. 3. The design procedure for the ZSS initial guess, based on the Kirchhoff-Love shell model, is described in Sect. 4. The numerical examples are given in Sect. 5, and the concluding remarks in Sect. 6.

EBTL method
We first provide, from [70], an overview of the EBTL method [67], including the EBZSS and IPBZSS concepts and the conversion between the two ZSS.
Let 0 ⊂ ℝ 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 ⊂ ℝ 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, is the displacement, is the virtual displacement, is the variation of the Green-Lagrange strain tensor, is the second Piola-Kirchhoff stress tensor, 0 is the mass density in the ZSS, is the body force per unit mass, and is the external stress vector applied on the subset t h of the boundary t . (1)

EBZSS
In the EBTL method the ZSS is defined with a set of positions 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, REF ∈ REF , all elements are connected by nodes, and we measure the displacement from that connected state. The implementation of the method is simple. The deformation gradient tensor is evaluated for each element: where 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 e 0 to integration-point counterpart of 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 , can be expressed as The Green-Lagrange strain tensor, where 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 IJ can be expressed with the components of the metric tensors. Thus, the inner product ∶ , and all the other quantities, in the weak form given by Eq. (5) can be evaluated without actually using the basis vectors I . This justifies using G IJ k as the integration-point counterpart of e 0 , with k = 1, … , n int , where n int is the number of integration points. Note that G IJ 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 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 IJ k , we solve a steady-state element-based problem (with = and = ): and the solution to that, in the form REF + , 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 A to be zero and constrain B to be in the direction The last constraint is C to be on the plane defined by the vector

Coordinate systems for the artery inner surface and wall
A geometrical relationship between the ZS and reference states, for a straight tube, was described in [70]. It was based on the shell model, where it is assumed that the inner-surface elements are extruded in the normal direction. Here we increase the scope of the method to general parametrization in the computational space. For each integration point of the computational space, we use a special shell-like coordinate system specific to that integration point. We will explain that coordinate system later in this section, after first explaining how the mesh is generated. In our notation here, will now imply REF , which is our "target" shape, and will imply 0 . We explain the method in the context of one element across the wall. Extending the method to multiple elements is straightforward.

Mesh generation
We again start with the artery inner surface, and build the wall in some fashion. Here we first build a T-spline inner-surface mesh. Then we expand that by an estimated thickness to obtain the outer-surface mesh. After that we modify the outer-surface mesh manually by moving the control points. When the thickness is larger than the radius of curvature, parts of the outer surface overlap. This might happen near the branches. Therefore, a simple extrusion does not work. Since the outer surface cannot be obtained from medical images, the design of the outer surface has to be based on other anatomical knowledge. This is the reason why currently meshes are generated manually rather than by an automated process. After defining the outer-surface mesh, which will have a control point corresponding to every control point on the inner surface, we add two control points for each pair. We use the four control points to form a cubic Bézier element across the wall. This is the way we obtain a T-spline volume mesh.

Inner-surface coordinates in the target state
The natural coordinates is used for the surface representation as explained in [70]. With the notation • indicating the inner surface, the basis vectors are where = 1, 2 . We define the unit normal vector as and the second fundamental form as and the curvature tensor is We also define two orthonormal vectors 1 and 2 , which are the principal directions of the curvature tensor: where ̂1 and ̂2 are the corresponding eigenvalues, and ̂1 ≥̂2 . The derivation and more details are in [70].

Inner-surface coordinates in the ZSS
Since the principal-curvature directions 1 and 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 unit normal vector in the thickness direction is We assume that the two directions 1 and 2 are also the principal-stretch directions. Then the ZSS basis vectors are calculated from where ̂1 and ̂2 are the principal stretches. Based on that we obtain the ZSS covariant basis vectors as where t = ⋅ . The derivation and more details are in [70].

Wall coordinates in the target state
In a target volume element, we again use the natural coordinate system. The position in the target state is ( ) , and the total differential of the position is Here we introduce a coordinate in the "normal" direction: where ̂ is the unit normal vector inside the wall. As we perform the integration ∫ d , the difference between the values we find as we reach the inner and outer surfaces will be our formal definition of the thickness h th . There are two options for defining ̂ .
The first option is that we find the closest point ⊥ on the inner surface from a given position : assuming a reasonable geometry for the purpose of calculating ̂ . For ‖ ( ) − ⊥ ‖ = 0 , we get Eq. (30). Figure 1 shows the normal definition. The second option is which naturally gives Eq. (30) at the surface. Figure 2 shows the normal definition. The other two components of the new coordinate system is represented by the vector ̂∈ ℝ n sd −1 , and the corresponding basis vectors are ̂ : For the normal-vector definition of Eq. (42), it simplifies to Remark 1 We note that even if ̂ = , ̂ and are not the same in general, because 3 is not perpendicular to 1 and 2 and will have an out-of-plane component.
The total differential of the position is expressed as

Remark 2
The first normal-vector option is closer to the shell theory. The downside is that we may not always have a reasonable geometry, for example when the radius of curvature is low compared to the thickness. The second option does not suffer from that problem and gives us the possibility of coming up with an ̂ design that would make the method better. However its quality depends on the mesh, such as the smoothness of the constant-3 surfaces.
At the kth integration point we define a shell-like coordinate system: where k is an offset defined in such a way that ( (̂, 0)) is on the surface. To find the value, we use the following expression:      By differentiating from Eq. (46) with respect to ̂ , we obtain the covariant basis vectors in the thickness direction: Here b is the second fundamental form: which can be calculated by using the natural coordinates: How to obtain this expression is described in Appendix A.
To obtain Î , we inner-product the right-hand sides of Eqs. (39) and (45) The derivative of the contravariant basis vector needed in Eq. (56) is derived in Appendix B.

Wall coordinates in the ZSS
Once we define the shell-like coordinate system for kth integration point, we can use the corresponding ZSS coordinate system as described in [70]. From that we compute the components of the metric tensor corresponding to the natural coordinates.
Here we redescribe the method by using the notation from the earlier parts of Sect. 3. 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. This thickness will be expressed as  More explicitly, we can write

Analytical expression and design for the ZSS initial guess
As explained in [70], the design parameters are the principal curvatures ̂0 1 and ̂0 2 , and the stretches ̂1 and ̂2 for each principal-curvature direction. However, they are not a set of independent values because there are constraints for the ZSS to give us the target shape for the given load, and we try to get close to that even for the initial guess. To that end, we use the force equilibrium in the normal direction by using a local analytical solution for each integration point. After that, we describe the design for the ZSS initial guess.

Analytical solution based on the Kirchhoff-Love shell model
We use the Kirchhoff-Love shell model with plane-stress condition to obtain the analytical solution. That is the generalized version of the solutions given in [71] for pressurized sphere and cylinder.
To deal with an arbitrary geometry, we assume the shape is given in terms of the principal curvatures. That is the reason we use only the equilibrium equation in the normal direction and focus on a small surface area . From that surface area, we extrude in the normal direction by h th to define a volume: = ∫ h th 0 ( )d . The force equilibrium in the normal direction is where and p are the Cauchy stress tensor and pressure load. Now we express this by using the shell-coordinate system. Using the plane-stress condition, we can write and the divergence of it is

what is left in the normal direction is
The small area can be expressed as That can be simplified as

Thus, Eq. (87) becomes
To calculate ̃ , we need the ZSS. For that, we substitute = into Eq. (38) and obtain the covariant basis vectors of the ZSS: and from 0 ( ) , yet to be calculated, we obtain the covariant basis vectors: From that, we obtain the contravariant basis vectors: Because we assume that the principal stretches are also in the principal-curvature directions, we obtain the principal stretch in direction as From that, we obtain To obtain 0 ( ) , we start from = 0 = 0 and integrate by using Eq. (69). This requires numerical integration, unless the constitutive law is very simple.

The design for the ZSS initial guess
As proposed in [69,70], 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 [72]. The wall thickness is about 8-10 % of the diameter at the target configuration. This means that ̂1h th is 0.16-0.20, and ̂2h th is nearly equal to zero but it could be negative. Patient-specific geometries are more complicated than that. To classify the local shapes, we consider the quadrants of the space formed by ̂1h th and ̂2h th . In the first quadrant, we have a balloon-like shape, which may be seen at a branch or an aneurysm. In the second and fourth quadrants, we have saddle points, which may be seen near branches. In the third quadrant, we have both curvatures negative, which may be seen also near branches. We note that because of the Kirchhoff-Love shell assumption, the following conditions (for = 1 or 2) are out of scope:

Analytical solution for constant stretch on the inner surface
We assume that the stretches ̂1 =̂2 = 1.05 and find the opening angles 1 and 2 in terms of ̂1h th and ̂2h th . We define an opening angle for each direction: where = 1, 2.

Remark 4
We note that the stretch value 1.05 gives 1 = 410 • for a straight tube with ̂1h th = 0.16 , and 1 = 310 • with ̂1h th = 0.20.  Figure 3 shows 1 , which is obtained iteratively and may not be unique. The figure also gives us 2 when we interchange ̂1h th and ̂2h th . We note that in many regions we have symmetry with respect to ̂1h th =̂2h th . We see large departure from the symmetry in regions where ̂1̂2 < 0 , which are the saddle points.
We use this map for generating the initial guess. When the geometry is out of scope, we set = 0.

Straight tube
The tube has ̂1h th = 0.20 . The meshes used in the computations are shown in Figs. 4 and 5. They are based on cubic B-splines, with 16, 8 and 1 elements in the circumferential, longitudinal and thickness directions. For the first mesh, 3 is in direction, and for the second mesh, 3 is skew to direction. As can be clearly seen in Fig. 5, the mesh is twisted. Figures 6 and 7 show, for the mesh in Fig. 4, the stretches from the ZSS initial guess and converged ZSS. Figures 8 and  9 do the same for the mesh in Fig. 5. From all these results we see that our initial guess is very good, even for the mesh where 3 is skew to direction.

Y-shaped tube
This geometry was motivated by branched arteries. The tube parts have ̂1h th = 0.16 , with constant thickness. Figures 10  and 11 show the shape and ̂1h th and ̂2h th . The geometry has two umbilical points, where ̂1 =̂2 , and three saddle-point regions. The mesh is shown in Fig. 12. Figure 13 shows the IPBZSS, using the EBZSS representation, from the ZSS initial guess and converged ZSS. The stretches are shown in Figs. 14 and 15. The initial guess has nearly uniform stretch on the inner surface. However, at the saddle point, the maximum stretch is higher than what we obtain from the ZSS initial guess. This affects the stretch on the outer surface too. There is less effect at the umbilical points. Figure 16 shows the stretch in ̂ direction. Again, there is less effect at the umbilical points.

Patient-specific aorta
The largest diameter is about 30 mm . Figures 17 and 18 show the shape and ̂1 and ̂2 . The inner-surface mesh is shown in Fig. 19. We use Laplace's equation to determine a smooth thickness distribution, setting the values at the boundaries to result in ̂1h th = 0.20 . We then generate the volume mesh as described in Sect. 3.1, which involves modification of the outer surface and consequently the thickness. The volume mesh is shown in Fig. 20. The measured thickness is displayed on the inner-surface mesh in Fig. 21, with an average value of ̂1h th = 0.18.
From the results of the Y-shaped tube we learned that the initial guess should have had, on the outer surface, less variation in the stretch in ̂ direction. Based on that, we improve the initial guess in parts of the region where we do not have both curvatures positive. For a straight tube with ̂1h th = 0.18 , on the outer surface, the stretch in ̂ direction is 0.80. We target that value in improving the initial guess.
We define four cases in deciding how to do the improvement, which we go through in the order given below:  Figure 25 shows the stretch in ̂ direction. Overall, the ZSS initial guess and converged ZSS are very similar. This indicates that reaching the ZSS design targets based on the analytical solution works well.

Concluding remarks
We have increased the scope and robustness of the method we introduced earlier for estimating the ZSS required in patient-specific arterial FSI computations, where the medical-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 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 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.
The method has two main components. 1. An iteration technique, which starts with a calculated ZSS initial guess, is used for computing the IPBZSS such that when a given pressure load is applied, the medical-image-based target shape is matched. 2. A design procedure, which is based on the Kirchhoff-Love shell model of the artery, is used for calculating the ZSS initial guess.  Fig. 19 Patient-specific aorta. Inner-surface mesh, made of cubic and quartic T-splines. Red circles represent the control points. The parts with the quartic T-splines, obtained by order elevation [73], are around the two extraordinary points, each connected to six edges We increased the scope and robustness of the method by introducing a new design procedure for the ZSS initial guess. The new design procedure has two features. (a) An IPB shell-like coordinate system, which increases the scope of the design to general parametrization in the computational space. (b) Analytical solution of the force equilibrium in the normal direction, based on the Kirchhoff-Love shell model, which places proper constraints on the design parameters. This increases the estimation accuracy, which in turn increases the robustness of the iterations and the convergence speed.
To show how the new design procedure for the ZSS initial guess performs, we first presented 3D test computations with a straight tube and a Y-shaped tube. After that we presented a 3D computation where the target geometry is coming from medical image of a human aorta, and we included the branches in the model. The computations show that the new design procedure for the ZSS initial guess is reaching the design targets well. We start with the transformation from the contravariant basis vectors to the covariant basis vectors: We take the derivative of both sides: and from that obtain From that and using g = ⋅ , we obtain Multiplying both sides with g , we obtain Thus,

Appendix C: Fung's model
The elastic-energy density function for the Fung's model is where D 1 and D 2 are the coefficients of the Fung's model. As a compressible-material model, we use the following expression: (117) = g − ⋅ .