A new approach to calculating fiber fields in 2D vessel cross sections using conformal maps

: An arterial vessel has three layers, namely, the intima, the media and the adventitia. Each of these layers is modeled to have two families of strain-sti ff ening collagen fibers that are transversely helical. In an unloaded configuration, these fibers are coiled up. In the case of a pressurized lumen, these fibers stretch and start to resist further outward expansion. As the fibers elongate, they sti ff en, a ff ecting the mechanical response. Having a mathematical model of vessel expansion is crucial in cardiovascular applications such as predicting stenosis and simulating hemodynamics. Thus, to study the mechanics of the vessel wall under loading, it is important to calculate the fiber configurations in the unloaded configuration. The aim of this paper is to introduce a new technique of using conformal maps to numerically calculate the fiber field in a general arterial cross-section. The technique relies on finding a rational approximation of the conformal map. First, points on the physical cross section are mapped to points on a reference annulus using a rational approximation of the forward conformal map. Next, we find the angular unit vectors at the mapped points, and finally a rational approximation of the inverse conformal map is used to map the angular unit vectors back to vectors on the physical cross section. We have used MATLAB software packages to achieve these goals.


Introduction
Cardiovascular diseases, leading to heart attacks (myocardial infarction) and stroke, are among the prime causes of mortality in modern society.Often increased stress levels, lack of exercise and many other factors are associated with cardiovascular disorders.Myocardial infarction is one of the most common causes of fatality and occurs when blood flow to the heart muscles is blocked.It has been the subject of many investigations over the past few decades [19], [16], [23], [20], [11].
One of the most important factors behind myocardial infarction is the blockage of arteries (stenosis).As coronary arteries undergo blockage, the heart experiences more pressure to maintain the blood flow, resulting in increased stress on the heart muscles, usually leading to heart failure.Thus, it is of prime importance to provide a mathematical model that can capture the dynamics of stenosis, and a crucial component of these models is an accurate mechanical model of the artery [4].
For the mechanical model, we need a strain energy density function from which we can calculate the lowest energy deformation; from the strain energy function, the stress response of a hyperelastic material can be derived (see [5], Chapter 6).The mathematical form of the strain-energy density function for nearly incompressible arterial tissue consists of several terms: a neo-Hookean term and a contribution from collagen fibers [3].In general, the strain energy density function takes the form W = W NH (I 1 ) + W aniso (I 1 , I 4 , I * 4 ) (1.1) where W NH (I 1 ) = µ(I 1 − 3) is a neo-Hookean function of I 1 = Trace(C) with µ being a layer-specific material parameter, C = F T F is the right Cauchy-Green tensor that registers the deformation from the reference to the deformed domain, and F is the deformation gradient.A commonly used term for W aniso (I 1 , I 4 , I * 4 ) is The term W aniso (I 1 , I 4 , I * 4 ) registers the anisotropy of the model [7,17] due to the presence of two families of collagen fibers with normalized lengths I 4 and I * 4 (see Section 1.1 below for more details).The constants ρ, η and β refer to the material constants of each of the layers.The other prime factor to consider while taking into account the dynamics of stenosis is the stiffness of the collagen fibers.In Eq (1.2), βρ determines the stiffness of a layer.A higher value of βρ implies a higher degree of stiffness compared to a lower value of βρ.The degree of stiffness could dictate the rate of lumen occlusion due to intimal growth [10].
In a healthy coronary artery, the intima consists of a single layer of smooth muscle cells (SMCs) and no collagen fibers.Over time there is an influx of SMCs from the media into the intima.A gradual thickening of the intima occurs when the SMCs undergo cell division.The proliferation of the SMCs in the intima and accumulation of lipids result in the formation of lesions.The authors in [8] observed localization of collagen fibers around the lesions.
Although there is no histological evidence of transverse helical collagen fibers in the intima, studies like [3,7] have assumed a two-family fiber model for the strain energy functional for the intima.
Over the past few decades, some intensive investigations have been carried out to simulate and provide models for the distribution of fibers in the layers of an arterial vessel.When it comes to the myocardium, fibers also play an important role, and there have been many studies that attempt to numerically construct fiber fields within the myocardial wall.In [15], the authors introduce a Laplace-Dirichlet Rule-Based (LDRB) algorithm to assign fiber orientations to heart models.In this study, the authors use linear interpolation for calculating fiber orientations within the myocardium.Deterministic techniques like interpolation are very common when modeling fibers in the heart.In [24], the authors solve a Poisson equation subjected to the Dirichlet conditions that the fiber orientations in the heart vary gradually from −70 • to 80 • from the outer to the inner wall.The authors in [21] study a post-mortem human heart.They reconstruct a continuous fiber field from experimental measurements of the actual fiber direction at a set of discrete points within each slice.In [1], the authors perform a diffusion tensor imaging (DTI) of seven carotid plaques.They observe that in advanced stage carotid plaques, the predominant fiber orientation is in the circumferential direction, with the elevation angle being zero or close to zero.
Fibers also play an important role in governing the dynamics and mechanics of an arterial wall.Microscopically, fibers are oriented helically around the lumen and are mostly aligned with wall boundaries with a small amount of dispersion [14].Therefore, in studies of 2D annular cross sections, fiber fields are assigned to the circumferential direction by simply computing the unit tangential vector.However, in reality, arterial cross sections are not perfectly annular.In such cases, the fiber orientations are sometimes modeled using interpolation techniques similar to the ones discussed above for myocardial tissue.For example, the authors in [3] calculated the fiber orientations as weighted averages of the tangent vectors on the Jordan boundaries of an arterial layer.
In addition to several interpolation models for computing and simulating fiber orientations, there are several studies on probabilistic approaches for computing fiber field orientations.In [4], structural continuum frameworks are developed, taking into account the dispersion of the collagen fiber orientations.A probability density function of the fiber orientation in the reference frame is used in [6] for the construction of the strain energy function.Another study [18] used maximum-likelihood estimation to determine the angle of the fibers along with a model to quantify collagen fiber dispersion.Similarly, in [9], a density function describing the distribution of the collagen fibers is incorporated into the strain energy function, and it evolves as the fibers orient themselves in the direction of the maximum Cauchy stress.The authors in [2] establish a relationship between the fiber end points and their arc lengths, which is estimated from image data.The mechanical response of collagen fibers, subjected to stress, is determined by a random variable approximated by a beta distribution.There is also open source software, FEBio (https://febio.org/),that simulates fiber orientations using probability distributions.In FEBio, a random distribution of fiber orientation is assigned to the points on the domain, and the average strain energy density function is calculated (https://help.febio.org/FebioUser/FEBioum 3-4-Section-4.3.html).The type of distribution function can be chosen from a range of distribution functions, which provides one with options of selecting spherical, ellipsoidal, circular, elliptical, von Mises and π-periodic von Mises distributions according to the geometry under consideration (https://help.febio.org/FebioUser/FEBioum 3-4-Subsection-4.3.3.html).
It is of significant importance to model the orientation of the fiber fields and predict how they change in order to describe accurately the response of the arterial walls due to loading.Computational techniques to explicitly compute fiber fields in general arterial cross sections using non-probabilistic approaches are rarely documented.In this study, we introduce conformal maps as an alternative to interpolation techniques when it comes to modeling fiber fields in 2D vessel cross sections.

Microstructural considerations
We consider an infinitely long cylindrical artery that has its cross section aligned with the X-Y plane, and every cross section is mechanically and geometrically identical for every Z.Each of the three layers of an artery (intima, media and adventitia) has transverse helical collagen fibers associated to it.These collagen fibers offer resistance to the expansion of the arteries.In the reference configuration, they remain in their lowest energy state.As soon as the arteries start to dilate under pressure, the collagen fibers start to uncoil, offering resistance to further dilation of the artery.As the collagen fibers stiffen gradually, it takes more and more energy to dilate the vessel.The collagen fibers have an inclination angle.Inclination angle is the angle made by the transversely helical collagen fibers with respect to the X-Y plane.To incorporate the effect of the collagen fibers on deformation, we introduce a variable I 4 for each of the layers, intima, media and adventitia.I 4 is the normalized length of the collagen fibers.As the vessel is loaded, I 4 increases.However, in reality, there are two families of fibers.One family is oriented at an angle ϕ, with respect to the axial direction, while the other is oriented at an angle π − ϕ.We introduce I 4 and I * 4 to represent fiber lengths for each family: T are the unit fiber vectors at the point (X, Y) in the reference stress-free/unloaded domain corresponding to each of the families of the fibers.In our work, since we ignore the residual stresses, the stress-free configuration is the same as the unloaded configuration.From the definitions of b and b * , we see that I 4 = I * 4 .Hence, W aniso in Eq (1. 2) only depends on I 1 and I 4 and takes the form The inclination angle, ϕ, is different for each layer, so the length of the fiber, projected onto the X-Y plane, is also different for each layer.In the reference stress-free/unloaded configuration, I 4 = 1 since C = 1, and b is a unit vector.Fb gives us the deformation of the unit fiber vector in the deformed configuration.

Motivation behind conformal maps
The function b(X, Y) for any physical cross section is not easy to calculate for a general point (X, Y).In the previous section, we explored different strategies for computing b(X, Y).Since the calculation is done numerically, many algorithms prove to be computationally intensive.The main idea here is to come up with a new approach using numerical conformal maps to calculate b(X, Y).
We introduce a reference annulus with inner radius ρ and outer radius 1, with 0 < ρ < 1.Given a physical cross section, the idea is to construct a conformal map that maps the physical cross section to the reference annulus.There are a few reasons behind choosing conformal maps.First, the fiber field in the reference annulus is very easy to compute.For any point (U, V) on the annulus, the unit fiber direction vector can be derived easily by calculating the vector − sin(θ) cos(θ) T , where θ = [atan2(V, U)]; atan2 is the four-quadrant arctangent.The other advantage of a conformal map is that it preserves angles locally.The fiber fields in a general cross section should satisfy some basic properties.First, fiber vectors at the boundaries between the intima, media and adventitia should be tangential: We assume that the fibers from one layer do not penetrate the adjacent layers.Second, the field should vary smoothly across the layers of the vessel.A smooth fiber field is essential for the strain-energy function to be well defined.General maps that are not conformal will not preserve angles: They can map tangential vectors in the annulus to non-tangential vectors in the physical cross section.However, for conformal maps, vectors tangent to boundaries remain tangential.Thus, the aim is to create a conformal map from the reference annulus to the physical cross section in Figure 1.Since the fiber field in the reference annulus can be computed easily, we can find the fiber field in the physical cross section simply by multiplying the unit fiber vector δω with the derivative of the inverse conformal map: where δω = − sin θ +i cos θ is the unit fiber vector on the reference annulus at the point ω, δz is the fiber vector in the physical cross section, and g is the inverse conformal map from the reference annulus to the physical cross section.The projected vector is proportional to [Re(δz) Im(δz)] T .This vector is inclined at an angle ϕ.Therefore, b ∝ Re(δz) Im(δz) Re 2 (δz) + Im 2 (δz) tan ϕ T .Therefore, in order for b to have unit length, we can write where L = (1 + tan 2 ϕ)(Re 2 (δz) + Im 2 (δz)), and ϕ is the inclination angle of the fibers in each of the layers.In the study [7], the authors inferred the inclination angles (along with other material constants) by fitting to stress-strain data.They showed that in human coronary arteries, ϕ ≈ 60.Conformal mapping of fiber fields between the reference annulus and the physical cross section.f is the forward conformal map from the physical cross section to the reference annulus, and g is the inverse conformal map from the reference annulus to the physical cross section.We map z to ω using f , find the unit fiber vector δω in the reference annulus and map it back to the fiber vector δz in the physical cross section using the derivative of the inverse conformal map g.

Numerical construction of the conformal map
This section will focus on the mathematical framework and the construction of the conformal map numerically using MATLAB.

The mathematical framework
Let's consider simply connected domains as a first step to understanding the method.The idea [22] is to construct the unique conformal map f : Ω → D, with Ω being a smooth simply connected domain bounded by a Jordan curve P and D being the open unit disk, such that f (0) = 0 and f ′ (0) > 0.Then, h(z) = log( f (z)/z) is a non-zero analytic function on Ω that is continuous on Ω and can be written as h(z) = u(z) + iv(z), where u(z) and v(z) are real harmonic functions in Ω and continuous on Ω.Any point z p = re iθ on the Jordan boundary P is mapped to f (z p ) = e iψ , where ψ ∈ R, with | f (z p )| = 1.Then, u(z p ) + iv(z p ) = i(ψ − θ) − log r, where r, θ ∈ R for z p ∈ P.
Writing f (z) in terms of h(z), we get where u(z) satisfies the Laplace equation with a Dirichlet boundary condition: In the above equation, v(z) is a harmonic conjugate of u(z) satisfying v(0) = 0. Solving Eq (2.2) and finding its harmonic conjugate v(z) gives the algebraic form of the conformal map f (z) according to Eq (2.1).The solution to Eq (2.2) is approximated by the real part of a complex analytic function: for some n, where c k = a k + ib k , a k , b k ∈ R, and b 0 = 0, Therefore, Since z ∈ P and u(z) = − log |z|, we plug in distinct points sampled from P into − log |z|, forming the left-hand side of Eq (2.4).For the right-hand side of Eq (2.4), we get a weighted linear combination of a k and b k for each of the points z ∈ P. To set up a least squares problem, let's sample 8n points on the Jordan curve P, z 0 , z 1 , . . ., z 8n−1 .Then, from Eq (2.4), we get  The above matrix equation can be written in a compact form as Ax ≈ d.Solving the least squares problem in Eq (2.5), we get the coefficients a 0 , a 1 , . . ., a n and b 1 , b 2 , . . ., b n .Plugging the coefficients into Eq (2.4), we get u(z), and since v(z) is the harmonic conjugate of u(z), it takes the form (2.6) Once we get u(z) and v(z), they can be substituted in Eq (2.1) to get the conformal mapping f (z).For a smooth doubly connected domain, there are two Jordan boundaries P 1 and P 2 , and we employ a similar strategy to the one above.However, there are differences in the formulation of the problem [22].The strategy for a smooth doubly connected domain is to find a conformal map f (z) that maps the domain to a circular annulus with inner radius ρ such that 0 < ρ < 1, where ρ is known as the conformal modulus, and the outer radius is 1.If P 1 and P 2 are the outer and the inner Jordan boundaries, respectively, of the smooth doubly connected domain, then f satisfies f (P 1 ) = S and f (P 2 ) = ρS , where S is the unit circle.Note that ρ is an unknown that must be found.The approximation of u(z) for a smooth doubly connected domain takes the form where c k = a k + ib k , for some a k , b k ∈ R and b 0 = 0. Since v(z) is the harmonic conjugate of u(z), we get (2.9) Thus, the Dirichlet boundary conditions for u(z) will be (2.10) Similar to finding u(z) for a smooth simply connected domain, we solve the Laplace equation ∆u = 0 subjected to Dirichlet boundary conditions in Eq (2.10) to find u(z) and its harmonic conjugate v(z).Since there are two Dirichlet boundary conditions for u(z), we get two matrix equations.To set up a least squares problem, we sample 8n points from each of the Jordan boundaries.From P 1 , we sample the points z 0,1 , z 1,1 , z 2,1 , . . ., z 8n−1,1 .Similarly, for the Jordan boundary P 2 , we sample the points z 0,2 , z 1,2 , z 2,2 , . . ., z 8n−1,2 .Now, since log ρ is unknown, the conformal map can be found by solving the least squares system Solving the least squares problem in Eq (2.11), we get the coefficients as well as the conformal modulus.From the least squares solution, we can construct u(z) and v(z) to get the conformal map f (z) according to Eq (2.1).We implement the conformal map numerically using the open source package chebfun (https://www.chebfun.org/,http://www.chebfun.org/examples/)and conformal2 (http://www.chebfun.org/examples/complex/ConformalMapping2.html).
The Fourier series approximations of Jordan boundaries P 1 and P 2 take the parametric forms r 1 (θ) and r 2 (θ), respectively.Multiplying the parametric forms by e iθ gives us the sampled points on the Jordan boundaries in the forms r 1 (θ) cos(θ) + ir 1 (θ) sin(θ) for P 1 and r 2 (θ) cos(θ) + ir 2 (θ) sin(θ) for P 2 .The function chebfun samples 8n points from each of r 1 (θ)e iθ and r 2 (θ)e iθ to set up the least squares problem in Eq (2.11).The function conformal2 solves the least squares problem to find the series approximation of f (z), and then the AAA algorithm takes the series expansion form of f (z) and converts it to a barycentric rational approximation [13].The reason for using the AAA algorithm is that the rational approximation uses a smaller number of coefficients than a polynomial representation of f (z) [13].
The AAA algorithm samples m, with m ≤ 16n, distinct points, also known as support points, from the Jordan boundaries of the physical cross section.Let's call these points z 1 , z 2 , . . ., z m .These points are then mapped to their corresponding images, f 1 = f (z 1 ), f 2 = f (z 2 ), . . ., f m = f (z m ), on the boundaries of the reference annulus using the forward map f (z) in Eq (2.1), using a set of non-zero complex weights w 1 , w 2 , . . ., w m .The function conformal2 returns a function handle f, with the values of z j , f j and w j , for j = 1, 2, 3, . . ., m.Then, the barycentric rational representation of f (z) takes the form (2.12) (see Theorem 2.1 [13]).According to [13], for reasonable choice of z j , the construction of the barycentric rational approximation becomes well conditioned.This improves the numerical stability of the conformal map.
For the inverse conformal map g, the AAA algorithm samples k distinct points, Z 1 , Z 2 , . . ., Z k , from the boundaries of the reference annulus, with k ≤ 16n.For these points, the corresponding points F j on the Jordan boundaries of the physical cross section are calculated, satisfying g(Z j ) = F j for j = 1, 2, . . ., k, using a set of non-zero complex weights W 1 , W 2 , . . ., W k .The code conformal2
returns another function handle for the inverse conformal map g(z), written as finv in the codes in the GitHub link (provided towards the end of this paper), containing the values of Z j , F j and W j for j = 1, 2, 3, . . ., k, from which the barycentric rational approximation of the inverse conformal map can be constructed in the form g(z) = N(z)/D(z).Then, the derivative of g(z) can be written as (2.13) Then, the unit fiber vector in the physical cross section can be found using Eq (1.5).The code conformal2 also returns the value of the conformal modulus ρ.

Results and discussion
Given a general arterial cross section, we find the forward and the inverse conformal map through barycentric rational approximation of f and g.Thus, for each of the layers, we find the fiber field on the corresponding reference annulus and map the fiber field back to the physical cross sections by applying Eqs (1.5) and (1.6).Following the mapping procedure as shown in Figure 1, the resulting fiber fields are computed in Figure 2.  Conformal mapping of fiber fields between the reference annulus and the physical cross section for each of the three layers of an artery.Here, we show mappings for two different cross sectional geometries.ρ is the conformal modulus.The values of m for the conformal map for the intima in (d) to (a), media in (d) to (b) and adventitia in (d) to (c) are 29, 30 and 26, respectively.The values of m for the conformal map for the intima in (h) to (e), media in (h) to (f) and adventitia in (h) to (g) are 34, 24 and 25, respectively., respectively.On the media-adventitia interface, we define the inclination angle to be 67 • , and on the intima-media interface, we define the inclination angle to be 60.3 • .
We also applied our technique to Figure 2(a) of [12], which is a cross section of an adult right coronary artery, and generated the fiber fields in each of the three layers.

Conclusions
The computation of fiber fields is of utmost importance in an accurate model for arterial mechanics.We have shown that, by using conformal maps, we can achieve that goal efficiently.To calculate the fiber fields using conformal maps, we only need a representation of the Jordan boundaries of the three layers.Provided with the Fourier series, the MATLAB function conformal2 then calculates the conformal modulus ρ, the forward function handle f and the inverse function handle finv, which are necessary to create the barycentric rational approximations of the forward and the inverse maps.Using the barycentric rational approximations, we then create the fiber fields for each of the layers in the physical cross section.
However, there are limitations in using the method of conformal maps.Our model does not take into account fiber dispersion as probabilistic approaches do.However, that can be readily fixed by introducing appropriate distribution functions to the fiber orientation in the reference or the physical cross section.Another limitation is that the code conformal2 diverges for Jordan boundaries that have points with large second derivatives (i.e., where the curvature is very high).This commonly results in the crowding phenomenon [22] (see the GitHub link below for an example).
As discussed in the Introduction, collagen fibers in the myocardium or an artery are mainly oriented in a "circumferential direction"; however, the notion of a circumferential direction is not well defined in the case of a general cross section geometry.The term only makes sense in the case of an annular cross section.
Conformal maps make the concept of a circumferential direction mathematically precise.We say that a fiber vector is oriented in a circumferential direction in a general cross section if and only if the fiber when mapped to the reference annulus using the forward conformal map is the angular unit vector.
Our method works for any general cross section geometry (except for the one limitation mentioned above).Given a cross section geometry, once the map is constructed, the fiber directions can be found out readily.The lengths of the fibers projected onto the physical cross section depend on the angles of inclination.This results in fibers of different projected lengths in different layers.Also, computational techniques for explicitly finding the fiber field for general cross section geometries using non-probabilistic approaches are rarely documented.This paper fills the gap and provides an example of an important biomechanical application of conformal mapping.
Figure 1.Conformal mapping of fiber fields between the reference annulus and the physical cross section.f is the forward conformal map from the physical cross section to the reference annulus, and g is the inverse conformal map from the reference annulus to the physical cross section.We map z to ω using f , find the unit fiber vector δω in the reference annulus and map it back to the fiber vector δz in the physical cross section using the derivative of the inverse conformal map g.

Figure 2 .
Figure 2.Conformal mapping of fiber fields between the reference annulus and the physical cross section for each of the three layers of an artery.Here, we show mappings for two different cross sectional geometries.ρ is the conformal modulus.The values of m for the conformal map for the intima in (d) to (a), media in (d) to (b) and adventitia in (d) to (c) are 29, 30 and 26, respectively.The values of m for the conformal map for the intima in (h) to (e), media in (h) to (f) and adventitia in (h) to (g) are 34, 24 and 25, respectively.

Figure 2 (
Figure 2 (a)-(c) represents the fiber field in the reference annulus for each of the three layers.(a) represents the reference annulus for the intima of the physical cross section in (d).Similarly, (b) and (c) represent the reference annuli for the media and the adventitia, respectively, corresponding to the physical cross sections in (d).(e)-(h) represent the reference annuli for the intima, media and adventitia and the physical cross section, respectively, for a different physical cross section geometry.The conformal modulus ρ is different for each of the reference annuli.A 3D visualization of the fibers is shown in Figure 3.

Figure 3 .
Figure 3. 3D view of the fiber fields for the three layers.Region colored red is the intima, region colored green is the media, and the region colored brown is the adventitia.The angles of inclination with respect to the plane for the intima, media and the adventitia are 60.3 • , 20.61 • and 67 • , respectively.On the media-adventitia interface, we define the inclination angle to be 67 • , and on the intima-media interface, we define the inclination angle to be 60.3 • .

Figure 4 .
Figure 4. Fiber fields generated by the method of conformal maps for the cross section in Figure2(a) of[12].