**1. Introduction**

Any perturbations in the blood stream, which is a biologically active fluid flowing in a channel with biologically active walls, inevitably lead to the activation of the body's defense systems and/or damage to the walls of the flow channel. Therefore, the main hypothesis is that the blood flow in the central parts of the circulatory system (heart and grea<sup>t</sup> vessels) is carried out without the formation of separation and stagnant zones, that is, it is a potential flow. In such a flow, by definition, any types of interaction are minimized both in the flow core and on the walls of the flow channel.

The cellular composition of the walls [1,2], their biomechanical properties [1], the distribution of velocities and shear stresses [3,4], and metabolism [5,6] were studied. Many model studies of the aorta have been carried out [7]. However, by present days there are no proper quantitative criterions that allow one to formalize the degree of aortic duct pathological remodeling. In the present study these parameters have been obtained by considering analytical solutions for the velocity vector of swirling blood flow in the heart and grea<sup>t</sup> vessels.

### **2. Materials and Methods**

In previously published works [8–10], it has been shown that the dynamic geometry of the heart and grea<sup>t</sup> vessels corresponds with a high degree of accuracy to the direction

**Citation:** Talygin, E.; Gorodkov, A.; Tibua, T.; Bockeria, L. Quantitative Criteria for the Degree of Pathological Remodeling of the Aortic Duct. *Mathematics* **2022**, *10*, 4773. https://doi.org/10.3390/ math10244773

Academic Editor: Mauro Malvè

Received: 30 October 2022 Accepted: 10 December 2022 Published: 15 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of streamlines of swirling flows described by exact solutions of the Navier-Stokes and continuity equations for a class of self-organizing tornado-like flows of a viscous fluid. These solutions were obtained in 1986 by G.I. Kiknadze and Yu.K. Krasnov [11] and are generally expressed by relation (1).

$$\begin{cases} u\_r = \mathbb{C}\_0(t)r + \frac{\mathbb{C}\_1}{r} \\ u\_z = -2\mathbb{C}\_0(t)z + \mathbb{C}\_2(t) \\ u\_\phi(r, t) = \frac{1}{r} \ast \left( A\_1 + A\_2 \Gamma\left(1 + \frac{\mathbb{C}\_1}{2v'}, a(0) \ast r^2\right) \right), \end{cases} \tag{1}$$

Here *<sup>C</sup>*0(*t*) and *<sup>C</sup>*2(*t*) are time-dependent functions, *A*1, *A*2, and *C*1 are constants, and *α*(*t*) is a function that has the following form:

$$\kappa(t) = \frac{e^{-2\int\_0^t \mathbb{C}\_0(\tau\_1)d\tau\_1}}{B\_1 - A \int\_0^t e^{-2\int\_0^t \mathbb{C}\_0(\tau\_2)d\tau\_2}d\tau\_1}$$

Here *B*1—is the arbitrary constant.

In the steady state, the swirling flow under consideration can be described by the relations for the Burgers vortex [12]:

$$\begin{cases} \quad \mu\_r = -\mathbb{C}\_0 \ast r\\ \quad \mu\_z = 2\mathbb{C}\_0 \ast z\\ \mu\_\theta = \frac{\Gamma\_0}{2\pi r} \ast \left(1 - e^{-\frac{\xi\_0^2 \ast r}{2r}}\right) \end{cases} \tag{2.1}$$

However, it is known that the blood flow is roughly unsteady and occurs in a pulsating regime. Previously, it has been shown [13,14] that the geometry of the flow channel during the entire cardiac cycle corresponds to the direction of the streamlines described by these solutions. Therefore, we used these solutions as quasi-stationary, if only *<sup>C</sup>*0(*t*) and <sup>Γ</sup>0(*t*) depend on time. In this case, relation (2.1) may be transformed as follows:

$$\begin{cases} \begin{array}{c} \boldsymbol{\mu\_{r}} = -\mathsf{C}\_{0}(t) \ast \boldsymbol{r} \\ \boldsymbol{\mu\_{z}} = 2\mathsf{C}\_{0}(t) \ast \boldsymbol{z} \\ \boldsymbol{\mu\_{\boldsymbol{\theta}}} = \frac{\Gamma\_{0}(t)}{2\pi r} \ast \left(1 - e^{-\frac{\mathsf{C}\_{0}^{2}(t) \ast r}{2\boldsymbol{v}}}\right) \end{array} \tag{2.2}$$

Here, only the azimuthal velocity component depends on the viscosity; therefore, the modulus of this component decreases with the evolution of the flow.

The geometry of the flow channel must correspond to the direction of the streamlines of the swirling flow inside. Expression (2.2) includes two functions of time—*C*0(*t*) and <sup>Γ</sup>0(*t*). To use these solutions as quasi-stationary in a non-stationary pulsating flow, the cardiac cycle was considered as a sequence of discrete states, in each of which the flow is described by relations (2.1). In this case, the instantaneous values of *<sup>C</sup>*0(*t*) and <sup>Γ</sup>0(*t*) are described by the measured geometric characteristics of the flow channel.

The instantaneous position of the streamlines of the considered flow in the longitudinalradial projection is described by the following relation:

$$zr^2 = \beta(\mathbf{t})\tag{3}$$

*β*(*t*)—a time-dependent function, the instantaneous value of which is determined by the choice of values *C*0 and Γ0 at a given time.

Streamlines in a fixed position of the flow channel of the heart and aorta in the axialradial projection are described by the following expression:

$$\varphi = \frac{\Gamma\_0}{2\pi \mathcal{C}\_0} \ast \left( \frac{1}{r} + \sqrt{\frac{\pi \mathcal{C}\_0}{2\upsilon}} \ast \int\_0^{\sqrt{\frac{\mathcal{C}\_0}{2\upsilon}}r} e^{-t^2} dt - \frac{1}{r} \ast e^{-\frac{\mathcal{C}\_0 r^2}{2\upsilon}} \right) + \varphi\_0 \tag{4}$$

As has been shown already, the aorta retains the shape of a converging canal throughout the entire cardiac cycle [13]. According to our assumptions, the pulse oscillations of the aortic wall are much less than the geometric transformations that occur in the process of remodeling. Therefore, the aortic flow channel was normally considered as a tube with almost constant geometric characteristics. As a result, fluctuations in the values of the characteristics *C*0 and Γ0 were considered negligible.

Since the streamlines are an evolution of the hyperbolic helix (4), continuous flow along such streamlines is possible only in a channel that has the form of a lower order hyperbolic right-handed helix (with fewer turns per aortic length). To approximate the shape of the aortic flow channel, a flat hyperbolic spiral of the following form was used:

$$r = \left(\frac{a}{\varphi}\right)^{\text{power}} + bias$$

Here (*r*, *φ*) are the radial and longitudinal coordinates, and (*a*, *power*, *bias*) are the parameters by which the approximation is carried out.

As a result, the obtained quantitative characteristics of the shape of the aorta are composed of two components, which are the values taken by the constants *C*0 and Γ0 and the parameters describing the helical properties of the aortic flow channel.

Let us consider 2 sections of the aortic flow channel outlying at a distance *z*1 and *z*2 respectively from the plane of the aortic valve. The corresponding aortic radiuses on these sections are *r*1 and *r*2. For *z*1 < *z*2 it follows that *r*1 > *r*2 due to the convergence of the flow channel. Since the specific spatial location of the point of origin of the cylindrical coordinate system in which the parameters of the aortic flow channel were calculated is unclear, a fictitious point of origin was introduced into consideration, which is located at a distance z0 from the plane of the aortic valve inside the cavity of the left ventricle. This point can change its position during the entire cardiac cycle, however, within the framework of this work, it was assumed to be immobile. It can be considered as the point of origin of the swirling blood flow. Then relation (3) for the sections under consideration can be written as follows:

$$(z\_0 + z\_1) \* r\_1^2 = (z\_0 + z\_2) \* r\_2^2 \Leftrightarrow z\_0 = \frac{z\_2 r\_2^2 - z\_1 r\_1^2}{r\_1^2 - r\_2^2} \tag{5}$$

From (2.2) it can be derived that the contribution of the azimuthal component in the total swirling flow velocity vector decreases with the distance from the origin point. Thus, if a certain section of the aortic flow channel *S*1 is closer to the point of origin than the section *S*2, the corresponding value of the modulus of the azimuth component can be expressed as follows: *<sup>u</sup>φ*,<sup>1</sup> > *<sup>u</sup>φ*,2. Accordingly, the minimum value of *<sup>u</sup>φ* occurs near the end of the aorta. For simplicity it can be stated as *<sup>u</sup>φ* = 0. Then the total velocity vector through the section *S*2 at the end of the aorta will be written as follows:

$$\mu\_2 = \sqrt{\mu\_z^2 + \mu\_r^2} = \mathbb{C}\_0 \sqrt{4(z\_0 + z\_2)^2 + r\_2^2} \Leftrightarrow \mathbb{C}\_0 = \frac{\mu\_2}{\sqrt{4(z\_0 + z\_2)^2 + r\_2^2}} \tag{6}$$

Considering the section *S*1, the derived expression for the linear velocity of the swirling flow through this section can be written as follows:

$$\begin{split} \mu\_1 &= \sqrt{\mathfrak{u}\_z^2 + \mathfrak{u}\_r^2 + \mathfrak{u}\_\rho^2} \Leftrightarrow \mathfrak{u}\_\rho^2 = \mathfrak{u}\_1^2 - \mathfrak{u}\_z^2 - \mathfrak{u}\_r^2 \Leftrightarrow \left( \frac{\Gamma\_0}{2\pi r} \ast \left( 1 - e^{-\frac{\mathfrak{L}\_0^2 \ast \mathfrak{r}\_1}{2\nu}} \right) \right)^2 \\ &= \mathfrak{u}\_1^2 - 4\mathbb{C}\_0^2 \left( z\_0 + z\_1 \right)^2 - \mathbb{C}\_0^2 r\_1^2 \end{split}$$

Expanding the brackets and simplifying the expression above, we get:

$$\Gamma\_0 = \frac{\sqrt{4\pi r\_1^2 u\_1^2 - 4\pi r\_1^2 \mathcal{C}\_0^2 \left(4(z\_0 + z\_1)^2 + r\_1^2\right)}}{\left(1 - e^{-\frac{\mathcal{C}\_0^2 \ast r\_1}{2\nu}}\right)}\tag{7}$$

For the correct reconstruction of a parametric spatial curve approximating the shape of the aortic flow channel, it is necessary to reconstruct the central line of the aorta first.

The central line is a spatial curve drawn between two planes cutting through an extended cavity in such a way that the distance from this line to the boundaries of the cavity is the maximum.

To construct the central line, 2 planes were used that cut the aortic flow channel perpendicular to the direction of the swirling blood flow. For all studied aortas, one of these planes is the plane that includes the aortic valve. The algorithm for choosing another plane depends on the type of aorta. For a relatively healthy aorta, the second plane should coincide with the section of the aorta in a bifurcation zone; for an aorta with a pathological disorder of the vascular bed, the plane was chosen approximately at the level of 2/3 of the aorta length, counting from the aortic valve. This is due to a large distortion of the geometry of the aortic flow channel in the abdominal region, which is associated with serious errors in the reconstruction of the required central line. Then, points were fixed on the selected planes in the central region of the cutting planes (one for each plane). These points (labeled as *p*1 and *p*2) were used to construct the required center line.

To determine the central line, it is necessary to determine the trajectory *C* = *<sup>C</sup>*(*s*) connecting the selected points *p*1 and *p*2. This trajectory should ensure the minimization of the value of the following functional:

$$I\_{centerline}(\mathbb{C}) = \int\_{\mathbb{C}^{-1}(p\_1)}^{\mathbb{C}^{-1}(p\_2)} F(\mathbb{C}(s)) ds \tag{8}$$

In the written expression, *<sup>F</sup>*(*x*) is a certain scalar field, the value of which at points lying closer to the center of the cavity is less than at points far from the center. The simplest example is a function whose value at a point is inversely proportional to the distance from this point to the boundary of the cavity. Such a function can be represented by the following relation:

$$DT(\mathbf{x}) = \min\_{y \in d\Omega} (|\mathbf{x}\_\prime y|) \tag{9}$$

In this relation, |(*<sup>x</sup>*, *y*)| is the Euclidean distance from the point x to the point y, *d*Ω is the boundary of the cavity Ω, corresponding to the radius of the channel at this point. Choosing the scalar field *<sup>F</sup>*(*x*) = *DT*−<sup>1</sup>(*x*) defined by expression (8), the central lines will need to be located on the medial axis of the given cavity Ω. The medial axis of the cavity is defined as the set of ball centers inscribed in a given cavity, where the size of the inscribed ball will be the maximum only if it in turn is not inscribed in another such ball. The construction of such a medial axis and the central line associated with it is carried out using the Voronoi diagram [15]. In the three-dimensional case, the Voronoi diagram is a surface composed of convex polygons whose vertices are the centers of the maximum inscribed balls.

The construction of the Voronoi diagram was performed by Delaunay triangulation [16], which was accompanied by the removal of polygons that partially fell out of the given cavity. This method allowed us to reformulate the problem described by expression (6) in the form of an eikonal equation, a non-linear partial differential equation:

$$\left| \nabla T(\mathbf{x}) = F(\mathbf{x}) \right| \tag{10}$$

The boundary condition for the written equation is *<sup>T</sup>*(*p*1) = 0. Such an equation can be solved using the fast sweep method. As soon as the solution of equation (8) was obtained for the entire Voronoi diagram, the backpropagation method was used to construct a trajectory from point *p*1 to point *p*2 in the direction of the maximum decrease in the value of the scalar field (i.e., in the direction of the medial axis of the cavity).

Central lines have been constructed using The Vascular Modelling Toolkit (vmtk).

All required data engineering and math computations has been done using standard python packages—pandas, numpy, scipy, scikit-learn and vmtk.
