Wall Shear Stress Topological Skeleton Analysis in Cardiovascular Flows: Methods and Applications

: A marked interest has recently emerged regarding the analysis of the wall shear stress (WSS) vector ﬁeld topological skeleton in cardiovascular ﬂows. Based on dynamical system theory, the WSS topological skeleton is composed of ﬁxed points, i.e., focal points where WSS locally vanishes, and unstable/stable manifolds, consisting of contraction/expansion regions linking ﬁxed points. Such an interest arises from its ability to reﬂect the presence of near-wall hemodynamic features associated with the onset and progression of vascular diseases. Over the years, Lagrangian-based and Eulerian-based post-processing techniques have been proposed aiming at identifying the topological skeleton features of the WSS. Here, the theoretical and methodological bases supporting the Lagrangian- and Eulerian-based methods currently used in the literature are reported and discussed, highlighting their application to cardiovascular ﬂows. The ﬁnal aim is to promote the use of WSS topological skeleton analysis in hemodynamic applications and to encourage its application in future mechanobiology studies in order to increase the chance of elucidating the mechanistic links between blood ﬂow disturbances, vascular disease, and clinical observations.


Introduction
Recent advances in medical imaging, modeling, and computational fluid dynamics (CFD) have allowed the modeling of local hemodynamics in realistic, personalized cardiovascular models, fostering understanding of the association between local hemodynamics and the initiation and progression of vascular disease, and in a wider perspective, contributing to the translation of computational methods in real-world clinical settings to complement clinical information.
It has long been recognized that hemodynamic factors regulate several aspects of vascular pathophysiology [1,2]. Wall shear stress (WSS), the frictional force per unit area exerted by streaming blood on the endothelium, has been identified as a major biomechanical factor involved in vascular homeostasis. In fact, WSS is sensed through several vascular mechanosensors and biochemical machineries that regulate the expression of genes coding for extra-and intra-cellular proteins, playing a relevant role in the development, growth, remodeling, and maintenance of the vascular system [3,4]. In this scenario, a multitude of WSS-based descriptors of the near-wall hemodynamics has been proposed over the years to provide potential indicators of flow disturbances associated with aggravating biological events. In particular, regions at the luminal surface presenting with low [5] and oscillatory [6] WSS have been identified as localizing factors of vascular disease [3,6]. However, the complex hemodynamic milieu the endothelium is exposed to can be only partially characterized by low and oscillatory WSS [7,8], as confirmed by a large body of literature reporting poor-to-moderate (and sometimes, contradictory) associations between low and oscillatory WSS with respect to vascular disease, e.g., [7][8][9][10][11][12][13]. This indicates a limited current understanding of the mechanistic link between WSS and vascular disease that hampers the use of WSS not only as a biomarker of vascular disease but also as a predictor of its progression within a clinical context [14].
Stimulated by the need to improve the understanding of the link between altered hemodynamics and clinical observations, the topological skeleton of the WSS vector field at the luminal surface of an artery is receiving increasing interest [15][16][17][18][19][20]. Based on dynamical system theory, the WSS topological skeleton is composed of a collection of fixed points, i.e., focal points where WSS locally vanishes, and unstable/stable manifolds, consisting of contraction/expansion regions linking fixed points. Such an interest arises from the ability of WSS topological skeleton features to reflect cardiovascular flow features like flow stagnation, separation and recirculation that are known to be promoting factors for vascular disease [2,17]. Very recent studies have demonstrated that the WSS topological skeleton governs the near-wall biochemical transport in arteries [15,16,18,21], which plays a fundamental role in, e.g., the initiation of atherosclerosis and thrombogenesis [22,23]. In addition, evidence of a direct association between WSS topological skeleton features and markers of vascular diseases from real-world clinical data have recently emerged [20,24].
In the present study, we report and discuss the theoretical background of Lagrangianand Eulerian-based methods currently applied to the analysis of the WSS topological skeleton. Based on the recent promising findings highlighting a link between WSS topological skeleton features and markers of vascular disease [17][18][19][20][21]24], the aim of this study is to encourage the application of WSS topological skeleton analysis to cardiovascular flows as an ad hoc instrument that is potentially able to further elucidate the mechanistic link between WSS and vascular pathophysiology.

Topological Skeleton of a Vector Field
Topological features of a vector field have been largely studied in the context of dynamical systems theory. A dynamical system is defined as a set of n differential equations: .
x(t) = u(x, t); where t ∈ R + is the time, x 0 ∈ R n the initial position at time point t 0 , i.e., x 0 = x(t 0 ), and u(x, t) the velocity field. Given the initial condition x 0 ∈ R n , a unique solution of Equation (1) exists, called trajectory, given by: u(x(s), s) ds.
Associated with the dynamical system defined in Equation (1), the so-called flow map can be defined as follows: providing the expression of all the system trajectories at time t. In general, the topological skeleton of the vector field u is recognized to provide the organizing structures of the system itself.
In steady-state conditions (i.e., when vector field u(x, t) in Equation (1) does not explicitly depend on time), the topological skeleton of a vector field consists of a collection of fixed points ( Figure 1A) and the associated stable and unstable manifolds connecting them ( Figure 1B). A fixed point (or critical point) is a point x f p ∈ R n where the vector field locally vanishes. The nature of fixed points can be stable or unstable. A stable fixed point is characterized by a sink configuration, and it attracts the nearby trajectories, while an unstable fixed point is characterized by a source configuration, and it repels the nearby trajectories ( Figure 1A). A fixed point can be classified as a saddle point, node, or focus ( Figure 1A): (1) a saddle point is a point attracting and repelling nearby trajectories along different directions (i.e., where the streamlines of the vector field intersect themselves); (2) a stable/unstable node is characterized by converging/diverging streamlines; (3) a focus is characterized by spiraling trajectories, and it can be attracting or repelling.
Technically, the exact location of fixed points in a domain of interest can be identified by computing the Poincaré index [25], a topological invariant index quantifying how many times a vector field rotates in the neighborhood of a point. For the sake of simplicity, we consider the dynamical system in Equation (1) under steady-state conditions and lying in a 2D space, i.e., ( ) = ( ), ( ) , with ∈ ℝ . An explanatory example of how to calculate the Poincaré index can then be provided. Let ∈ ℝ be an isolated fixed point of with a neighborhood N such that there are no other fixed points in N than , and let γ be a closed curve inscribing N. Then, the Poincaré index ℐ(γ, ) of the curve γ relative to is the number of the positive field rotations while traveling along γ in a positive direction: where is the vector field rotation angle. The Poincaré index is equal to −1 at saddle point locations ( Figure 1A), 1 at node or focus locations ( Figure 1A), and 0 otherwise. The algorithm for computing the Poincaré index for a 3D vector field defined on unstructured triangle meshes is extensively described elsewhere [19].
The Poincaré index allows identifying fixed point locations, but it does not provide information about the fixed points nature. Therefore, a criterion to distinguish between a node or a focus and between the attractive or repelling nature of a fixed point is needed. In light of this, the vector field around the fixed point can be expressed by linearization as: where J is the Jacobian matrix of . The classification of fixed points can be thus performed by computing the eigenvalues of the Jacobian matrix J, as summarized in Table 1. In detail, two real eigenvalues with different signs identify a saddle point. Two real eigenvalues with the same sign identify nodes characterized as attracting or repelling (i.e., stable or unstable, respectively) according to their sign (negative or positive, respectively). Complex conjugate eigenvalues identify a stable or unstable focus according to the sign of the real part (negative or positive, respectively). A fixed point can be classified as a saddle point, node, or focus ( Figure 1A): (1) a saddle point is a point attracting and repelling nearby trajectories along different directions (i.e., where the streamlines of the vector field intersect themselves); (2) a stable/unstable node is characterized by converging/diverging streamlines; (3) a focus is characterized by spiraling trajectories, and it can be attracting or repelling.
Technically, the exact location of fixed points in a domain of interest can be identified by computing the Poincaré index [25], a topological invariant index quantifying how many times a vector field rotates in the neighborhood of a point. For the sake of simplicity, we consider the dynamical system in Equation (1) under steady-state conditions and lying in a 2D space, i.e., u(x) = (X(x), Y(x)), with x ∈ R 2 . An explanatory example of how to calculate the Poincaré index can then be provided. Let x f p ∈ R 2 be an isolated fixed point of u with a neighborhood N such that there are no other fixed points in N than x f p , and let γ be a closed curve inscribing N. Then, the Poincaré index I(γ, u) of the curve γ relative to u is the number of the positive field rotations while traveling along γ in a positive direction: where θ is the vector field rotation angle. The Poincaré index is equal to −1 at saddle point locations ( Figure 1A), 1 at node or focus locations ( Figure 1A), and 0 otherwise. The algorithm for computing the Poincaré index for a 3D vector field defined on unstructured triangle meshes is extensively described elsewhere [19]. The Poincaré index allows identifying fixed point locations, but it does not provide information about the fixed points nature. Therefore, a criterion to distinguish between a node or a focus and between the attractive or repelling nature of a fixed point is needed. In light of this, the vector field u around the fixed point x f p can be expressed by linearization as: where J is the Jacobian matrix of u. The classification of fixed points can be thus performed by computing the eigenvalues of the Jacobian matrix J, as summarized in Table 1. In detail, two real eigenvalues with different signs identify a saddle point. Two real eigenvalues with the same sign identify nodes characterized as attracting or repelling (i.e., stable or unstable, respectively) according to their sign (negative or positive, respectively). Complex conjugate eigenvalues identify a stable or unstable focus according to the sign of the real part (negative or positive, respectively).
The stable and unstable manifolds (or critical lines) associated with a fixed point x f p are composed of all initial conditions x 0 ∈ R n such that the trajectories initiated at these points x 0 approach the fixed point x f p asymptotically. By construction, stable and unstable manifolds act as separatrices of the vector field, portioning regions of different behavior and dynamics. In detail, an unstable manifold attracts nearby trajectories, as opposed to the stable manifold, which repels nearby trajectories ( Figure 1B). In mathematical terms, an unstable manifold W u associated with the generic fixed point x f p is defined as follows: while a stable manifold W s can be expressed as: In general, two different perspectives have been proposed to identify manifolds of a vector field, namely the Lagrangian and Eulerian perspectives. The Lagrangian perspective considers individual particles, tracking their motion along their paths as they are advected by the flow field. By contrast, the Eulerian perspective considers the properties of the vector field under analysis at each fixed location in space and time. In the following sections, a brief theoretical background is reported for a better understanding of the theory supporting the Lagrangian and Eulerian approaches for the analysis of vector field topology, with particular emphasis on their application to cardiovascular flows.

Lagrangian Coherent Structures
When the vector field u(x, t) in Equation (1) is time-dependent, solutions can be complex and chaotic, making the interpretation of the topological skeleton made of W u , W s and x f p difficult. The need to robustly define intrinsic structures governing fluid/mass transport under unsteady-state conditions has led to the development of the concept of coherent structures (CS). Technically, CS are defined as emergent patterns, influencing the transport of tracers/mass in time-dependent flows [26]. In this context, Lagrangian Coherent Structures (LCS) are coherent structures identified by applying methods based on a Lagrangian approach. The theoretical bases of LCS lie in methods of nonlinear dynamics, chaos theory, and fluid dynamics.
From a mathematical perspective and in relation to fluid mechanics, LCS can be defined as material surfaces in the flow field that are dominant in attracting or repelling neighboring fluid elements over a defined time interval [27,28]. These material surfaces are able to localize where the flow field experiences the largest and the smallest stretching [29]. In detail, material surfaces in the flow field attracting trajectories more strongly than any other nearby material surface are referred to as attracting LCS. Oppositely, material surfaces repelling trajectories more strongly than any other nearby material surface are referred to as repelling LCS.
The detection and visualization of LCS is usually performed by applying two different Lagrangian-based approaches, namely (1) Lagrangian particle tracking and (2) the computation of the finite-time Lyapunov exponent (FTLE). Both approaches are based on particle path information derived from the post-processing of velocity data obtained by CFD simu-Mathematics 2021, 9, 720 5 of 21 lation or by in vivo (e.g., phase contrast magnetic resonance imaging (MRI)) and in vitro (e.g., particle image velocimetry) measurements. The workflow of the Lagrangian-based approaches to visualize LCS is sketched in Figure 2.
Mathematics 2021, 9, x FOR PEER REVIEW 5 of 22 particle path information derived from the post-processing of velocity data obtained by CFD simulation or by in vivo (e.g., phase contrast magnetic resonance imaging (MRI)) and in vitro (e.g., particle image velocimetry) measurements. The workflow of the Lagrangianbased approaches to visualize LCS is sketched in Figure 2. The Lagrangian particle tracking is performed by seeding the domain of interest with tracer particles and by visualizing their motion ( Figure 2). The aim of this approach is to reveal coherent features revealing how the flow under analysis is organized. From a mathematical perspective, the position of a tracer particle is governed by the differential equation reported in Equation (1). To obtain the position of such a particle at a desired time , Equation (1) is numerically integrated from to . The direct integration of tracer particles allows for an in-depth understanding of how tracers are transported through the domain of interest. In detail, attracting LCS will be generally distinguishable, since tracer particles are attracted to and along these surfaces ( Figure 3). Analogously, repelling LCS will be distinguishable from the advection of tracer particles by reversing time (Figure 3). Attracting LCS are traced out with forward time integration of particles, while repelling LCS are traced out with backward time integration of particles. The Lagrangian particle tracking is performed by seeding the domain of interest with tracer particles and by visualizing their motion ( Figure 2). The aim of this approach is to reveal coherent features revealing how the flow under analysis is organized. From a mathematical perspective, the position of a tracer particle is governed by the differential equation reported in Equation (1). To obtain the position of such a particle at a desired time t, Equation (1) is numerically integrated from t 0 to t. The direct integration of tracer particles allows for an in-depth understanding of how tracers are transported through the domain of interest. In detail, attracting LCS will be generally distinguishable, since tracer particles are attracted to and along these surfaces ( Figure 3). Analogously, repelling LCS will be distinguishable from the advection of tracer particles by reversing time (Figure 3). Attracting LCS are traced out with forward time integration of particles, while repelling LCS are traced out with backward time integration of particles.  Lagrangian particle tracking represents a Lagrangian-based technique aiming at overcoming issues related to standard approaches used for topological skeleton extraction of vector fields with unsteady-state conditions. However, the resulting tracer particle motion complexity could obscure the interpretation of the vector field topology. For this motivation, the second approach consists of the computation of the FTLE ( Figure 2). Based on theory, a LCS can be defined as the material surface locally maximizing the FTLE [27,30], the Lyapunov exponent being a measure of the sensitivity to the initial position of a dynamical system. Technically, the finite time Lyapunov exponent ( , , ) [27,28,[31][32][33][34][35] is defined as: where ( , , ) is the maximum eigenvalue of the right Cauchy-Green strain tensor ( , , ): where ∇Φ ( ) denotes the transpose of the gradient of the flow map in Equation (3). From a physical point of view, ( , , ) in Equation (9) represents the material deformation of infinitesimal volume elements of fluid, and it is a symmetric and positive-definite matrix. Roughly speaking, the FTLE defined in Equation (8) measures the rate of separation of initially close vector field trajectories. Let be a small distance between two material points at time , as depicted in Figure 4 (Panel A). It can be demonstrated [26] that the separation after the time interval | − | satisfies the inequality: where equality holds if the initial distance is aligned with the eigenvector of ( , , ) associated with .
( , , ) = ∇Φ ( ) ∇Φ ( ), Lagrangian particle tracking represents a Lagrangian-based technique aiming at overcoming issues related to standard approaches used for topological skeleton extraction of vector fields with unsteady-state conditions. However, the resulting tracer particle motion complexity could obscure the interpretation of the vector field topology. For this motivation, the second approach consists of the computation of the FTLE ( Figure 2). Based on theory, a LCS can be defined as the material surface locally maximizing the FTLE [27,30], the Lyapunov exponent being a measure of the sensitivity to the initial position of a dynamical system. Technically, the finite time Lyapunov exponent σ(x 0 , t 0 , t) [27,28,[31][32][33][34][35] is defined as: where λ max (C(x 0 , t 0 , t)) is the maximum eigenvalue of the right Cauchy-Green strain tensor C(x 0 , t 0 , t): where ∇Φ t t 0 (x 0 ) T denotes the transpose of the gradient of the flow map in Equation (3). From a physical point of view, C(x 0 , t 0 , t) in Equation (9) represents the material deformation of infinitesimal volume elements of fluid, and it is a symmetric and positive-definite matrix. Roughly speaking, the FTLE σ defined in Equation (8) measures the rate of separation of initially close vector field trajectories. Let δ 0 be a small distance between two material points at time t 0 , as depicted in Figure 4 (Panel A). It can be demonstrated [26] that the separation δ t after the time interval |t − t 0 | satisfies the inequality: where equality holds if the initial distance δ 0 is aligned with the eigenvector of C(x 0 , t 0 , t) associated with λ max . The algorithm for LCS identification based on FTLE computation starts with the initialization of a cluster of massless elemental particles at time t 0 over the domain of interest ( Figure 2). Then, particles are numerically integrated by the field in Equation (1) from t 0 to t, and their trajectories are calculated. The flow map Φ t t 0 (Equation (3)) is obtained from the final position of each particle trajectory at time t in the domain, and subsequently its gradient ∇Φ t t 0 (x 0 ) can be computed. For a structured grid like the one shown in Figure 4 (panel B), ∇Φ t t 0 (x 0 ) can be calculated by finite differencing, e.g., using central differencing as follows: Mathematics 2021, 9, x FOR PEER REVIEW 7 of 22 The algorithm for LCS identification based on FTLE computation starts with the initialization of a cluster of massless elemental particles at time over the domain of interest ( Figure 2). Then, particles are numerically integrated by the field in Equation (1) from to , and their trajectories are calculated. The flow map Φ (Equation (3)) is obtained from the final position of each particle trajectory at time in the domain, and subsequently its gradient ∇Φ ( ) can be computed. For a structured grid like the one shown in Figure 4 (panel B), ∇Φ ( ) can be calculated by finite differencing, e.g., using central differencing as follows: Once the flow map gradient is obtained, the Cauchy-Green strain tensor ( , , ) can be computed according to Equation (9).
Finally, the maximum eigenvalue ( ( , , )) and the FTLE ( , , ) can be computed according to Equation (8) (Figure 2). The obtained ( , , ) value for each particle is assigned to the particle position at time . This procedure is repeated, varying the time (e.g., within the cardiac cycle in cardiovascular applications) and aiming at providing the time series of FTLE values and ultimately the time history of LCS movements ( Figure 2). Positive integration times reveal repelling LCS in the FTLE field, while negative integration times reveal attracting LCS in the FTLE field.
In general, the computation of the spatial variation of the FTLE field requires the vector field to be interpolated in both time and space, and high-order integration and interpolation schemes are needed to ensure accuracy of results. Furthermore, the mesh used to compute the FTLE distribution over the domain of interest usually needs to be more resolved than the computational mesh for a more robust detection of LCS.

LCS Application to Intravascular Flows
Lagrangian-based approaches have been largely applied to identifying LCS in intravascular flows. Indeed, Lagrangian particle tracking has been massively applied to explore the complexity of intravascular flows, e.g., to provide a measure of stasis in idealized computational bifurcation models [36], or to study vortices generation and their potential role in thrombogenesis in idealized aneurysm models [37,38]. Several studies have applied particle tracking to identify flow disturbances in, e.g., carotid bifurcation models, contributing to providing a deeper understanding of the hemodynamics-driven processes Once the flow map gradient is obtained, the Cauchy-Green strain tensor C(x 0 , t 0 , t) can be computed according to Equation (9).
Finally, the maximum eigenvalue λ max (C(x 0 , t 0 , t)) and the FTLE σ(x 0 , t 0 , t) can be computed according to Equation (8) (Figure 2). The obtained σ(x 0 , t 0 , t) value for each particle is assigned to the particle position at time t 0 . This procedure is repeated, varying the time t 0 (e.g., within the cardiac cycle in cardiovascular applications) and aiming at providing the time series of FTLE values and ultimately the time history of LCS movements ( Figure 2). Positive integration times reveal repelling LCS in the FTLE field, while negative integration times reveal attracting LCS in the FTLE field.
In general, the computation of the spatial variation of the FTLE field requires the vector field to be interpolated in both time and space, and high-order integration and interpolation schemes are needed to ensure accuracy of results. Furthermore, the mesh used to compute the FTLE distribution over the domain of interest usually needs to be more resolved than the computational mesh for a more robust detection of LCS.

LCS Application to Intravascular Flows
Lagrangian-based approaches have been largely applied to identifying LCS in intravascular flows. Indeed, Lagrangian particle tracking has been massively applied to explore the complexity of intravascular flows, e.g., to provide a measure of stasis in idealized computational bifurcation models [36], or to study vortices generation and their potential role in thrombogenesis in idealized aneurysm models [37,38]. Several studies have applied particle tracking to identify flow disturbances in, e.g., carotid bifurcation models, contributing to providing a deeper understanding of the hemodynamics-driven processes underlying atherosclerosis onset/progression [39][40][41][42]. Moreover, particle tracking has been used to study the hepatic perfusion in the Fontan circulation [43,44], identify the optimal left ventricular assist device cannula outflow configurations [45], obtain a deeper understanding of the dynamics of embolic particles within arteries [46], and detect peculiar intravascular helical flow patterns in the aorta from in vivo, 4D-flow MRI data [47,48].
Regarding the FTLE-based analysis of the flow field, its extension to intravascular flows is relatively recent, motivated by the fact that LCS are determined by blood flow structures associated to adverse vascular events including flow stagnation, separation, and recirculation. Among the main contributions, here we mention that Shadden and Taylor [32] used LCS to quantify the extent of flow stagnation to determine where flow separated and to understand how flow was partitioned to downstream vasculature in computational hemodynamic models of large vessels. LCS have been proposed as a powerful method to capture vortex transport in blood flow. In this regard, Arzani and Shadden [33] used LCS to characterize the hemodynamics in abdominal aortic aneurysm (AAA) models, suggesting that AAA intravascular flow topology is dictated by systolic vortex propagation through the abnormal vessel. Arzani et al. [49] computed FTLE fields and associated LCS to capture a large coherent vortex in AAA computational models. Furthermore, LCS have been applied to identify left ventricle (LV) blood flow features during heart filling. In detail, Gharib et al. [50] used LCS to demonstrate the existence of a link between the vortex ring formation inside the LV and the ejection fraction. Charonko et al. [51] quantified the vortex ring volume by computing LCS from in vivo LV phase contrast MRI data of healthy and diseased patients. Töger et al. [52] extracted LCS from in vivo LV phase contrast MRI data to measure the vortex ring volume during LV rapid filling. The identification of attracting and repelling LCS from LV Doppler-echocardiography data was adopted as a criterion to discriminate between healthy and diseased patients [53]. Other studies applied LCS to characterize the flow field through heart valves. In particular, LCS were extracted to delineate the boundaries of the outflow jet downstream of aortic valves and used as a measure of the severity of the valve's stenosis [54,55]. In a very recent study [56], FTLE-based LCS detection on computational hemodynamics models of aortic bicuspid and mechanical heart valves was used to study mass transport processes that might be related to valve disease. The analysis of the fluid dynamics in the neighborhood of blood clots was another effective application of LCS to hemodynamics [57]. In addition, FTLE-based analysis was adopted to highlight the hemodynamic impact of flow diverter stents in the treatment of intracranial aneurysms [58,59].
We refer the interested reader to reference [31] for a broader, detailed overview of Lagrangian methods used in post-processing of velocity data in cardiovascular flows.

LCS Application to Near-Wall Flow Features
Recently, in the study of cardiovascular flows, the concept of LCS has been extended to analyze the near-wall flow topology, i.e., the topology of the flow field close to the luminal surface of arteries. The rationale is in the well-established involvement of near-wall mass transport in most of the processes concurring to determine vascular pathophysiology [5]: in the near-wall region, blood flow regulates the local biotransport processes and imparts mechanical shear stress on the endothelium (i.e., the WSS), which in turn regulates important developmental, homeostatic, and adaptive mechanisms in arteries, as well as susceptibility to and progression of atherosclerosis [1].
Based on theory, it has been demonstrated [60] that the WSS vector field can be scaled to provide a first-order approximation for the near-wall blood flow velocity vector field as follows: where u π ∈ R 3 is the near-wall velocity, τ ∈ R 3 represents the WSS vector field, µ is the dynamic viscosity, and δn is the distance from the wall where the velocity is evaluated. By construction, the vector field in Equation (12) is defined on the luminal surface of the vessel, and it represents the near-wall velocity, as the velocity is zero on the surface itself due to the no-slip condition. The LCS underlying theory described in Section 3.1 can be extended to analyze the near-wall flow topology by using the expression of near-wall velocity u π (given by Equation (12)) in Equation (1). Such near-wall Lagrangian structures, computed from the WSS vector field, are referred to as WSS LCS [15].
Computationally, WSS LCS can be identified on the luminal surface of the vessel by numerically integrating a high number of luminal surface tracer particles, applying the procedure described in the first part of Section 3.1. In detail, attracting and repelling WSS LCS can be traced out with forward and backward time integration of surface tracer particles based on the near-wall blood flow velocity (Equation (12)), respectively.
The recent interest in WSS LCS from the cardiovascular fluid mechanics research community was driven by WSS LCS ability to highlight blood flow features associated with vascular disease initiation and progression, like flow stagnation, separation, recirculation, flow impingement, and the interaction of vortex structures with the vascular wall [19,61,62] These blood flow features have been classified as "aggravating flow events", as they trigger biological processes leading to the development or progression of vascular disease [2,17]. An example of attracting WSS LCS on the luminal surface of a patient-specific computational hemodynamic model of carotid bifurcation is presented in Figure 5. Details on the carotid bifurcation hemodynamic modeling are reported elsewhere [9,14,20,63]. In this specific case, luminal surface tracer particles ( Figure 5A) are numerically integrated in forward time. The resulting LCS is located at the carotid bulb, a region characterized by flow disturbances (slow, recirculating blood flow) promoting atherosclerosis [2,4]. In detail, the attracting WSS LCS provides the boundary at the luminal surface of the slow vortex structure formed inside the carotid bulb ( Figure 5C, where the recirculation region is highlighted visualizing the streamlines of the cycle-average velocity vector field).
Computationally, WSS LCS can be identified on the luminal surface of the vessel by numerically integrating a high number of luminal surface tracer particles, applying the procedure described in the first part of Section 3.1. In detail, attracting and repelling WSS LCS can be traced out with forward and backward time integration of surface tracer particles based on the near-wall blood flow velocity (Equation (12)), respectively.
The recent interest in WSS LCS from the cardiovascular fluid mechanics research community was driven by WSS LCS ability to highlight blood flow features associated with vascular disease initiation and progression, like flow stagnation, separation, recirculation, flow impingement, and the interaction of vortex structures with the vascular wall [19,61,62] These blood flow features have been classified as "aggravating flow events", as they trigger biological processes leading to the development or progression of vascular disease [2,17]. An example of attracting WSS LCS on the luminal surface of a patient-specific computational hemodynamic model of carotid bifurcation is presented in Figure 5. Details on the carotid bifurcation hemodynamic modeling are reported elsewhere [9,14,20,63]. In this specific case, luminal surface tracer particles ( Figure 5A) are numerically integrated in forward time. The resulting LCS is located at the carotid bulb, a region characterized by flow disturbances (slow, recirculating blood flow) promoting atherosclerosis [2,4]. In detail, the attracting WSS LCS provides the boundary at the luminal surface of the slow vortex structure formed inside the carotid bulb ( Figure 5C, where the recirculation region is highlighted visualizing the streamlines of the cycle-average velocity vector field). In addition, the shear forces exerted by the streaming blood flow in the near-wall region on the endothelium affect biotransport processes, i.e., the transport of biochemicals through the subendothelial layer [22]. Biotransport is of paramount importance in many cardiovascular processes, including the initiation of atherosclerosis and thrombogenesis [23]. In general, cardiovascular mass transport is investigated in silico by coupling the governing equations of motion, the Navier-Stokes equations, with the advection-diffusion equation, given by: In addition, the shear forces exerted by the streaming blood flow in the near-wall region on the endothelium affect biotransport processes, i.e., the transport of biochemicals through the subendothelial layer [22]. Biotransport is of paramount importance in many cardiovascular processes, including the initiation of atherosclerosis and thrombogenesis [23]. In general, cardiovascular mass transport is investigated in silico by coupling the governing equations of motion, the Navier-Stokes equations, with the advection-diffusion equation, given by: where C is a non-dimensional concentration of the species transported in the domain, u is the fluid velocity vector, and D is the mass diffusion coefficient. However, high computational costs are associated with the class of numerical simulations used to accurately solve the near-wall transport and blood-wall transfer [64,65], making this approach expensive in hemodynamics applications. To overcome this limitation, and based on the well-established role that WSS plays in conditioning the permeability of the endothelium and the near-wall mass transport process, recent studies [15,18] have brilliantly demonstrated that WSS LCS can be used as a template for near-wall mass transport. This allows reduction of the computational effort needed to solve the full transport problem, represented by Equation (13) [15]. In particular, it has been demonstrated that attracting WSS LCS attract biochemicals, leading to high near-wall concentration in their neighborhood, whereas repelling WSS LCS have been shown to act as near-wall transport barriers [15,17].
In the context of cardiovascular flows, it has been recently demonstrated that attracting/repelling WSS LCS on the luminal surface of an artery match the unstable/stable manifolds of the cycle-average WSS vector field [15,18], defined as: where τ is the instantaneous local WSS value and T is the time duration of the cardiac cycle. Technically, the first step in the topological analysis of cycle-average WSS at the luminal surface of a vessel is the identification of WSS fixed points. The exact position of WSS fixed points can be identified by computing, e.g., the Poincaré index, as explained in Section 2. Then, the cycle-average WSS field around a fixed point x f p , according to Equation (5), by linearization can be expressed as: where J is the Jacobian matrix of τ (see Equation (14)). The identified WSS fixed points can be classified according to their nature (i.e., saddle, node, or focus, Figure 1A) by analyzing the eigenvalues of the Jacobian matrix J of τ (Table 1), as described in Section 2. Note that the WSS vector field is embedded in a three-dimensional space even if it lies in a two-dimensional space (the luminal surface of a vessel). To perform a two-dimensional analysis, two strategies are possible. In the first strategy, a projection of the vector field into two orthogonal directions (hence, in a two-dimensional space) is needed. In the second one, avoiding the projection of the vector field (and thus reducing the computational steps), a three-dimensional analysis is performed, thus obtaining three eigenvalues of the Jacobian matrix, with one of them having a value close to zero. Then, the eigenvalue-based analysis for the WSS fixed points classification considers only the two eigenvalues different from zero. Saddle-type fixed points are of particular interest, since typically a stable or unstable manifold starts from a saddle point and vanishes into a source or sink, respectively, as depicted in Figure 1B. Saddle point locations (where the Poincaré index is −1 and the eigenvalues are real with different signs) are perturbed along the positive eigenvector of J in two opposite directions, obtaining two initial conditions [18,61]. Unstable manifolds can be traced out by numerically integrating τ from these initial conditions in forward time until trajectories reach a stable fixed point (sink configuration) or leave the domain. Similarly, stable manifolds are delineated by integrating τ in backward time starting from the perturbation of saddle point locations along the negative eigenvector of J until trajectories reach an unstable fixed point (source configuration) or leave the domain.
An example of unstable manifolds of cycle-average WSS on the luminal surface of a patient-specific computational hemodynamic model of carotid bifurcation is presented in Figure 6. Details on carotid bifurcation hemodynamic modeling are reported elsewhere [9,14,20,63]. WSS fixed points were identified by computing the Poincaré index ( Figure 6A), and subsequently, unstable manifolds were traced out by applying Runge-Kutta 4-5 numerical integration schemes ( Figure 6B). By visual inspection of Figure 6C, it can be appreciated that cycle-average WSS unstable manifolds co-localize with attracting WSS LCS, confirming the capability of the latter to identify critical lines of the WSS field. and subsequently, unstable manifolds were traced out by applying Runge-Kutta 4-5 numerical integration schemes ( Figure 6B). By visual inspection of Figure 6C, it can be appreciated that cycle-average WSS unstable manifolds co-localize with attracting WSS LCS, confirming the capability of the latter to identify critical lines of the WSS field. The analysis of cycle-average WSS fixed points and manifolds has been applied to analyze cardiovascular flows. Arzani et al. [18] used WSS LCS from stable and unstable manifolds of cycle-average WSS on patient-specific computational hemodynamics models of AAAs, carotid arteries, cerebral aneurysms, and coronary aneurysms to characterize near-wall flow topology and biochemical transport. Farghadan et al. [16] used WSS topology and magnitude analysis to predict surface concentration patterns in cardiovascular transport problems by computing WSS LCS from manifolds of cycle-average WSS in image-based coronary and carotid artery models. Mahmoudi et al. [21] studied the near-wall transport of some of the prominent biochemicals contributing to the initiation and progression of atherosclerosis in computational hemodynamic models of the coronary artery, highlighting the strength of cycle-average WSS LCS as a template for luminal surface concentration and flux patterns of biochemicals transported with blood.
Summarizing, the Lagrangian approach for identifying near-wall topological features is schematized in Figure 7, where the link between attracting/repelling WSS LCS with unstable/stable cycle-average WSS manifolds, respectively, is highlighted. In addition, Figure 7 presents a brief summary of the link between Lagrangian-based near-wall flow topology and mass transport. For a more in-depth analysis, the interested reader can refer to recent literature [15,17,18,21] where the link between WSS LCS, cycle-average WSS manifolds, and biochemical transport in cardiovascular flows is unambiguously documented. The analysis of cycle-average WSS fixed points and manifolds has been applied to analyze cardiovascular flows. Arzani et al. [18] used WSS LCS from stable and unstable manifolds of cycle-average WSS on patient-specific computational hemodynamics models of AAAs, carotid arteries, cerebral aneurysms, and coronary aneurysms to characterize nearwall flow topology and biochemical transport. Farghadan et al. [16] used WSS topology and magnitude analysis to predict surface concentration patterns in cardiovascular transport problems by computing WSS LCS from manifolds of cycle-average WSS in image-based coronary and carotid artery models. Mahmoudi et al. [21] studied the near-wall transport of some of the prominent biochemicals contributing to the initiation and progression of atherosclerosis in computational hemodynamic models of the coronary artery, highlighting the strength of cycle-average WSS LCS as a template for luminal surface concentration and flux patterns of biochemicals transported with blood.
Summarizing, the Lagrangian approach for identifying near-wall topological features is schematized in Figure 7, where the link between attracting/repelling WSS LCS with unstable/stable cycle-average WSS manifolds, respectively, is highlighted. In addition, Figure 7 presents a brief summary of the link between Lagrangian-based near-wall flow topology and mass transport. For a more in-depth analysis, the interested reader can refer to recent literature [15,17,18,21] where the link between WSS LCS, cycle-average WSS manifolds, and biochemical transport in cardiovascular flows is unambiguously documented.

Volume Contraction Theory
From a Eulerian perspective, the volume contraction theory provides a simple alternative way to analyze the behavior of a dynamical system. Contrarily to Lagrangian-based approaches, the Eulerian perspective considers vector field properties at each point in space and time. The here-presented volume contraction theory, based on fluid mechanics and differential geometry, is focused on the temporal change of an elemental volume (of fluid, for the case of interest) in the phase space of a dynamical system (fluid flow, for the case of interest). Let ( ) be an arbitrary volume in the phase space of the dynamical system defined in Equation (1). Let ( ) be a closed surface enclosing the volume ( ), i.e., such that ( ) = ( ). ( ) evolves during the time interval resulting in a contraction or expansion of the volume, as depicted in Figure 8. The rate of volume variation, which we will call volume contraction rate in the following, can be expressed as follows as a consequence of the application of the Gauss theorem: where is the vector field defined in Equation (1) and is the unit normal to the surface S (Figure 8). Shrinking the near-wall volume to a point, it can be written: Equation (17) clearly shows that in the limit as approaches zero, the local value of vector divergence is equal to its total flux per unit volume.

Volume Contraction Theory
From a Eulerian perspective, the volume contraction theory provides a simple alternative way to analyze the behavior of a dynamical system. Contrarily to Lagrangian-based approaches, the Eulerian perspective considers vector field properties at each point in space and time. The here-presented volume contraction theory, based on fluid mechanics and differential geometry, is focused on the temporal change of an elemental volume (of fluid, for the case of interest) in the phase space of a dynamical system (fluid flow, for the case of interest). Let V(t) be an arbitrary volume in the phase space of the dynamical system defined in Equation (1). Let S(t) be a closed surface enclosing the volume V(t), i.e., such that S(t) = δV(t). S(t) evolves during the time interval dt resulting in a contraction or expansion of the volume, as depicted in Figure 8. The rate of volume variation, which we will call volume contraction rate in the following, can be expressed as follows as a consequence of the application of the Gauss theorem: where u is the vector field defined in Equation (1) and n is the unit normal to the surface S ( Figure 8). Shrinking the near-wall volume V to a point, it can be written: Equation (17) clearly shows that in the limit as V approaches zero, the local value of vector u divergence is equal to its total flux per unit volume.
In general, in non-conservative dynamical systems, the volume of phase space is not preserved, as it can contract or expand. Thus, trajectories tend toward a lower-dimensional subset of phase space. From Equation (17), the volume contraction rate Λ(x,t) of a ndimensional system, representing the rate of separation of infinitesimal close trajectories, can be obtained as: where tr J(u) is the trace of the Jacobian matrix of vector field u and λ i are its eigenvalues. Physically, the Jacobian matrix describes how a small change at a starting point x 0 propagates to the final point of the flow map Φ t t 0 (x 0 ) of Equation (3). In this sense, Equation (18) represents the sum of the Lyapunov exponents of Equation (8). In general, in non-conservative dynamical systems, the volume of phase space is not preserved, as it can contract or expand. Thus, trajectories tend toward a lower-dimensional subset of phase space. From Equation (17), the volume contraction rate Λ(x,t) of a n-dimensional system, representing the rate of separation of infinitesimal close trajectories, can be obtained as: where ( ) is the trace of the Jacobian matrix of vector field and are its eigenvalues. Physically, the Jacobian matrix describes how a small change at a starting point propagates to the final point of the flow map Φ ( ) of Equation (3). In this sense, Equation (18) represents the sum of the Lyapunov exponents of Equation (8).

Eulerian-Based Approach for WSS Topological Skeleton Identification
It has been recently demonstrated [19] that the application of the volume contraction theory to cardiovascular flows allows the analysis of the WSS topological skeleton on the luminal surface of a vessel through the direct calculation of the WSS divergence. Briefly, considering the expression of the near-wall blood flow velocity vector given in Equation (12) and substituting it in Equation (17), it follows that: Based on Equation (19), the WSS divergence gives practical information about the WSS topological skeleton. Note that in general, the WSS vector field defined at the luminal surface of a vessel is not conservative, even in the case of incompressible flows.
Contextualizing the physical meaning of Equation (19) to the study of the phenomena at the interface between blood flow and vessel wall, it can be observed that as the divergence represents the volume density of the outward flux of a vector field from an infinitesimal volume around a given point: • a local positive value of the divergence of the WSS field at the luminal surface means that locally shear forces exert an expansion action on the endothelium; • a local negative value of the divergence of the WSS field at the luminal surface means that locally shear forces exert a contraction action on the endothelium.
In general, the application of the volume contraction theory to the analysis of a dynamical system faces one limitation in cases where the distance between two neighboring trajectories increases/decreases in spite of a negative/positive value of divergence, respectively. As WSS divergence depends by construction upon the algebraic summation of the magnitude of the single gradients of WSS vector components, in some cases, it might fail

Eulerian-Based Approach for WSS Topological Skeleton Identification
It has been recently demonstrated [19] that the application of the volume contraction theory to cardiovascular flows allows the analysis of the WSS topological skeleton on the luminal surface of a vessel through the direct calculation of the WSS divergence. Briefly, considering the expression of the near-wall blood flow velocity vector u π given in Equation (12) and substituting it in Equation (17), it follows that: Based on Equation (19), the WSS divergence gives practical information about the WSS topological skeleton. Note that in general, the WSS vector field defined at the luminal surface of a vessel is not conservative, even in the case of incompressible flows.
Contextualizing the physical meaning of Equation (19) to the study of the phenomena at the interface between blood flow and vessel wall, it can be observed that as the divergence represents the volume density of the outward flux of a vector field from an infinitesimal volume around a given point: • a local positive value of the divergence of the WSS field at the luminal surface means that locally shear forces exert an expansion action on the endothelium; • a local negative value of the divergence of the WSS field at the luminal surface means that locally shear forces exert a contraction action on the endothelium.
In general, the application of the volume contraction theory to the analysis of a dynamical system faces one limitation in cases where the distance between two neighboring trajectories increases/decreases in spite of a negative/positive value of divergence, respectively. As WSS divergence depends by construction upon the algebraic summation of the magnitude of the single gradients of WSS vector components, in some cases, it might fail to properly identify WSS expansion/contraction configuration patterns. In fact, these regions describe specific directional arrangements of the vectors, but both variations in magnitude and in direction are accounted for in the divergence. To overcome this limitation, which could markedly affect the application of the Eulerian-based approach to study WSS manifolds in cardiovascular flows, the use of the divergence of the normalized WSS vector field has been recently proposed [19]: where τ u is the WSS unit vector. Equation (20) can be used to encase the connections between fixed points, i.e., manifolds, identify basins of attraction, and subdivide the domain into different vector field behaviors. Then, in the light of the above and as depicted in Figure 9, luminal surface regions characterized by negative values of DIV W are referred to as contraction regions and approximate unstable manifolds. Similarly, regions at the luminal surface characterized by positive values of DIV W are referred to as expansion regions and approximate stable manifolds (Figure 9).
gions describe specific directional arrangements of the vectors, but both variations in mag-nitude and in direction are accounted for in the divergence. To overcome this limitation, which could markedly affect the application of the Eulerian-based approach to study WSS manifolds in cardiovascular flows, the use of the divergence of the normalized WSS vector field has been recently proposed [19]: where is the WSS unit vector. Equation (20) can be used to encase the connections between fixed points, i.e., manifolds, identify basins of attraction, and subdivide the domain into different vector field behaviors. Then, in the light of the above and as depicted in Figure 9, luminal surface regions characterized by negative values of are referred to as contraction regions and approximate unstable manifolds. Similarly, regions at the luminal surface characterized by positive values of are referred to as expansion regions and approximate stable manifolds (Figure 9). To complete the Eulerian-based WSS topological skeleton analysis, once WSS manifolds have been identified using , the WSS fixed point location can be carried out using the Poincaré index, as in the Lagrangian-based analysis (as described in Section 2). Then, the eigenvalues of the Jacobian matrix of the WSS vector field can be used to distinguish between a node or a focus and between the attractive or repelling nature of a fixed point, as described in Section 2 in general terms (Table 1) and in Section 3.3 for the specific case of a WSS vector field defined on the luminal surface of a vessel.
An example of Eulerian-based topological skeleton analysis of the cycle-average WSS field on the luminal surface of a patient-specific computational hemodynamic model of carotid bifurcation is presented in Figure 10A. Details on carotid bifurcation hemodynamic modeling are reported elsewhere [9,14,20,63]. WSS fixed points were identified and classified by computing the Poincaré index and eigenvalues of the Jacobian matrix, respectively, whereas contraction/expansion regions were identified by computing the divergence of the normalized cycle-average WSS vector field. By visual inspection of Figure  10B, it can be noted that cycle-average WSS contraction regions co-localize with cycleaverage WSS unstable manifolds, traced out by integrating cycle-average WSS starting from saddle point positions, thus confirming the capability of the contraction regions to encase WSS manifolds.
To complete the Eulerian-based WSS topological skeleton analysis, once WSS manifolds have been identified using DIV W , the WSS fixed point location can be carried out using the Poincaré index, as in the Lagrangian-based analysis (as described in Section 2). Then, the eigenvalues of the Jacobian matrix of the WSS vector field can be used to distinguish between a node or a focus and between the attractive or repelling nature of a fixed point, as described in Section 2 in general terms (Table 1) and in Section 3.3 for the specific case of a WSS vector field defined on the luminal surface of a vessel.
An example of Eulerian-based topological skeleton analysis of the cycle-average WSS field on the luminal surface of a patient-specific computational hemodynamic model of carotid bifurcation is presented in Figure 10A. Details on carotid bifurcation hemodynamic modeling are reported elsewhere [9,14,20,63]. WSS fixed points were identified and classified by computing the Poincaré index and eigenvalues of the Jacobian matrix, respectively, whereas contraction/expansion regions were identified by computing the divergence of the normalized cycle-average WSS vector field. By visual inspection of Figure 10B, it can be noted that cycle-average WSS contraction regions co-localize with cycle-average WSS unstable manifolds, traced out by integrating cycle-average WSS starting from saddle point positions, thus confirming the capability of the contraction regions to encase WSS manifolds.
The Eulerian-based approach to analyze the WSS topological skeleton can be easily implemented. It requires only single snapshots of the WSS vector field, and the postprocessing algorithms, based on a robust theory, are easily reproduced. This approach does not need the Lagrangian surface transport computation, as required for Lagrangian-based and integrated trajectory-based methods, thus reducing computational effort. Furthermore, it is characterized by modularity in the sense that the method can be applied only for the purpose of fixed point identification and/or classification or only for contraction/expansion region identification. The Eulerian-based approach to analyze the WSS topological skeleton can be easily implemented. It requires only single snapshots of the WSS vector field, and the post-processing algorithms, based on a robust theory, are easily reproduced. This approach does not need the Lagrangian surface transport computation, as required for Lagrangian-based and integrated trajectory-based methods, thus reducing computational effort. Furthermore, it is characterized by modularity in the sense that the method can be applied only for the purpose of fixed point identification and/or classification or only for contraction/expansion region identification.

Application of the Eulerian-Based Method for WSS Topological Skeleton Analysis to Cardiovascular Flows
The described Eulerian-based method to identify the WSS topological skeleton on the luminal surface of an artery can be easily applied (1) to cycle average WSS vectors (defined in Equation (14) in Section 3.3) and (2) to instantaneous WSS vectors along the cardiac cycle.
The cycle-average WSS topological skeleton highlights blood flow features associated with vascular disease development, and it is strongly related to arterial near-wall mass transport. In detail, on the one hand, contraction/expansion regions of cycle-average WSS vectors, because of their capability to encase unstable/stable cycle-average WSS manifolds, can be used to identify biochemical concentration patterns at the arterial luminal surface. On the other hand, the instantaneous WSS topological skeleton allows analyzing the unsteady nature of WSS fixed points and contraction/expansion regions. In detail, instantaneous WSS fixed points may have a potential impact on the endothelial cells (ECs) function. By definition, a WSS fixed point represents a focal point on the luminal surface of a vessel where WSS vanishes, and low WSS is a biomechanical factor involved in vascular dysfunction. In light of this, quantitative descriptors of WSS fixed points residence times along the cardiac cycle were proposed [17,19], aiming at characterizing their unsteady nature. In detail, a WSS fixed point residence time, that for each surface element measures the accumulated amount of time that WSS fixed points spend inside that element,weighted by the sum of the absolute values of the eigenvalues of the instantaneous WSS Jacobian matrix, was proposed elsewhere [17]. More recently, a different formulation for quantifying WSS fixed points was proposed [19] where the local residence time of WSS fixed points were weighted by the absolute value of the sum of the eigenvalues of the WSS

Application of the Eulerian-Based Method for WSS Topological Skeleton Analysis to Cardiovascular Flows
The described Eulerian-based method to identify the WSS topological skeleton on the luminal surface of an artery can be easily applied (1) to cycle average WSS vectors (defined in Equation (14) in Section 3.3) and (2) to instantaneous WSS vectors along the cardiac cycle.
The cycle-average WSS topological skeleton highlights blood flow features associated with vascular disease development, and it is strongly related to arterial near-wall mass transport. In detail, on the one hand, contraction/expansion regions of cycle-average WSS vectors, because of their capability to encase unstable/stable cycle-average WSS manifolds, can be used to identify biochemical concentration patterns at the arterial luminal surface. On the other hand, the instantaneous WSS topological skeleton allows analyzing the unsteady nature of WSS fixed points and contraction/expansion regions. In detail, instantaneous WSS fixed points may have a potential impact on the endothelial cells (ECs) function. By definition, a WSS fixed point represents a focal point on the luminal surface of a vessel where WSS vanishes, and low WSS is a biomechanical factor involved in vascular dysfunction. In light of this, quantitative descriptors of WSS fixed points residence times along the cardiac cycle were proposed [17,19], aiming at characterizing their unsteady nature. In detail, a WSS fixed point residence time, that for each surface element measures the accumulated amount of time that WSS fixed points spend inside that element, weighted by the sum of the absolute values of the eigenvalues of the instantaneous WSS Jacobian matrix, was proposed elsewhere [17]. More recently, a different formulation for quantifying WSS fixed points was proposed [19] where the local residence time of WSS fixed points were weighted by the absolute value of the sum of the eigenvalues of the WSS Jacobian matrix (i.e., according to Equation (18), the absolute value of the WSS divergence, representing the strength of the contraction/expansion of the WSS around the fixed point). In mathematical terms: where x f p is the location of a WSS fixed point at time t ∈ [0, T], T is the cardiac cycle duration, e is the generic triangular element of the superficial mesh of area A e , A is the average surface area of all triangular elements of the superficial mesh, I e is the indicator function (equal to 1 if x f p ∈ e, 0 otherwise) and (∇·τ) e is the instantaneous WSS divergence value, representing the local strength of the contraction/expansion of the WSS around the considered fixed point. Roughly speaking, Equation (21) allows quantifying the fraction of the cardiac cycle for which a generic triangle mesh surface element e on the vessel luminal surface hosted as a fixed point, weighting the residence time by the strength of the local contraction/expansion action. Furthermore, the strength and the nature of WSS contraction and expansion action on the ECs lining the luminal surface, as identified by WSS contraction/expansion regions, is expected to have biological consequences linked to vascular pathophysiology. In particular, the exposure to high variability of WSS contraction and expansion action may mechanically induce a recurring variation in EC stimulation along the cardiac cycle, with consequent widening cell-cell junctions and associated increased endothelium permeability and EC dysfunction and apoptosis [7,66]. The amount of variation in the WSS contraction/expansion action exerted at the luminal surface of a vessel along the cardiac cycle can be quantified using the quantity topological shear variation index (TSV I) [24]: Equation (22) allows localizing regions on the vessel luminal surface exposed to large variations in WSS contraction/expansion action exerted by the flowing blood along the cardiac cycle.
An example of the distribution of WSS fixed point weighted residence time (Equation (21)) and the topological shear variation index (Equation (22)) on the luminal surface of a patient-specific computational hemodynamic model of carotid bifurcation is presented in Figure 11. Details on the carotid bifurcation hemodynamic modeling are reported elsewhere [9,14,20,63]. From Figure 11, it emerged that the highest RT∇ x f p (e) values and highest variation in the contraction/expansion action exerted by the WSS along the cardiac cycle were mainly located at the carotid bulb and around the bifurcation apex. Interestingly, very recent studies highlighted a link between WSS contraction/expansion variability along the cardiac cycle and aggravating biological events at the arterial wall. In particular, De Nisco et al. [24] applied the Eulerian-based approach for the analysis of the WSS topological skeleton for personalized computational hemodynamic models of ascending thoracic aorta aneurysm (ATAA) and healthy aorta, reporting that: (1) the different spatiotemporal heterogeneity characterizing the ATAA and healthy hemodynamics markedly reflect on their WSS topological skeleton features; (2) a link emerged between the variability of the contraction/expansion action exerted by WSS on the endothelium (as quantified by the ) along the cardiac cycle and ATAA wall stiffness. Morbiducci et al. [20] demonstrated in a longitudinal study integrating clinical data with computational hemodynamics that WSS topological skeleton features quantified by the TSVI independently predicted long-term restenosis after carotid bifurcation endarterectomy. Interestingly, very recent studies highlighted a link between WSS contraction/expansion variability along the cardiac cycle and aggravating biological events at the arterial wall. In particular, De Nisco et al. [24] applied the Eulerian-based approach for the analysis of the WSS topological skeleton for personalized computational hemodynamic models of ascending thoracic aorta aneurysm (ATAA) and healthy aorta, reporting that: (1) the different spatiotemporal heterogeneity characterizing the ATAA and healthy hemodynamics markedly reflect on their WSS topological skeleton features; (2) a link emerged between the variability of the contraction/expansion action exerted by WSS on the endothelium (as quantified by the TSV I) along the cardiac cycle and ATAA wall stiffness. Morbiducci et al. [20] demonstrated in a longitudinal study integrating clinical data with computational hemodynamics that WSS topological skeleton features quantified by the TSVI independently predicted long-term restenosis after carotid bifurcation endarterectomy.

Future Directions
The translation into clinical settings of the WSS topological skeleton is hampered by several barriers that add up to those affecting in general the translation of computational hemodynamics and the derived knowledge, as discussed elsewhere [67]. Specifically, as a first step, the analysis of the topological skeleton needs to be distilled into intuitive, clinically relevant criteria. To this aim, only semi-quantitative results are obtained from the definition of fixed points and stable/unstable manifolds, consisting of contraction/expansion regions. However, quantitative results can be obtained by focusing on specific features by using ad-hoc topological skeleton descriptors, such as the fixed point weighted residence time RT∇ x f p (e) (Equation (21)) or the topological shear variation index, TSVI (Equation (22)). Then, the definition of clinically relevant criteria based on the WSS topological descriptors require cut-off values for an effective translation into the clinic. These cut-off values need to be accurately defined and tested in terms of performance including accuracy, sensitivity, specificity, and positive predictive value, among others. Therefore, the determination of cut-off values requires adequate statistical power, obtained usually through multiple prospective, randomized trials. Moreover, the endpoint to be predicted should be clearly defined, as different endpoints correspond to different cut-off values.
In the perspective of an effective translation into the clinic of quantitative topological skeleton features, in a previous study [20], we proved that exposure to high values of both descriptors RT∇ x f p (e) and TSVI was correlated with intima-media thickness (a marker of vascular disease) at 60 month follow-ups in carotid bifurcations after carotid endarterectomy. To determine the cut-off values of the descriptors, a pooled distribution for each descriptor was calculated from 46 models of healthy carotid bifurcation. The 80th percentile of those distributions was then used. This approach allowed definition of the cut-off values for abnormally high values of RT∇ x f p (e) or the TSVI.
It is evident that cut-off values are specific to the vascular region and to the predicted endpoint and therefore cannot be extrapolated to other conditions. In the future, the continuous improvements in imaging and data acquisition, the increasing availability of computational power, and the development of more and more efficient and robust methodologies for blood flow modeling are expected to accelerate the translation into the clinic of the analysis of the WSS topological skeleton. Our paper aims to give the methodological basis to tackle these future efforts.

Conclusions
The need for the identification of hemodynamic coherent structures in blood vessels is dictated by the so-called hemodynamic risk hypothesis, suggesting a major role of flow disturbances in vascular pathophysiology [2]. The action of fluid forces on the endothelial mechanosensors and biochemical machinery has been historically explained in terms of WSS [3,4]. However, only moderate (and sometimes contradictory) associations between vascular disease and WSS-based descriptors have emerged to date, motivating a more in-depth analysis of the fluid near-wall transport phenomena. In this sense, the capability of the WSS topological skeleton to capture features reflecting cardiovascular flow complexity [17][18][19][20] and having a direct link to adverse vascular biological events has recently attracted a strong research interest. In this regard, recent studies have demonstrated that the cycle-average WSS topological skeleton governs the near-wall biochemical transport in arteries [15,16,18], a process linked to, e.g., endothelium-mediated vasoregulation, thrombosis, and atherosclerosis [23]. Furthermore, evidence about the role of WSS topological skeleton features in vascular pathophysiology emerged from very recent studies suggesting a direct link between WSS topological skeleton features and, e.g., aortic wall stiffness [24] and late restenosis in endarterectomized carotid arteries [20].
Motivated by the need to characterize more precisely the WSS phenotype(s) linked to aggravating biological events, here we provided an overview of the theoretical and methodological basis for analyzing the WSS topological skeleton in cardiovascular flows. In detail, the present study is intended to: (1) promote the application of WSS topological skeleton analysis to cardiovascular flows, aiming at elucidating the role that peculiar WSS features play in vascular pathophysiology; (2) facilitate the reproducibility and comparability of results from WSS topological skeleton analyses among different studies; (3) confirm its potential as a tool for increasing the chance of elucidating the mechanistic link between flow disturbance and clinical outcomes when applied to real-world clinical data.
Here, both WSS topological skeleton Lagrangian-and Eulerian-based methods currently adopted in the literature are presented. Lagrangian-based approaches start from the processing of Eulerian data, which represent the typical outputs of current in vivo (e.g., phase contrast MRI), in vitro (e.g., particle image velocimetry), and computational methods used for the investigation of cardiovascular flows. On the one hand, Lagrangian approaches are particularly useful for revealing the global organization of the vector field and characterizing its evolution over time, making the relevant features easy to detect by visual inspection, as they offer effective three-dimensional (or even four-dimensional, i.e., including time) visualizations. On the other hand, Lagrangian techniques rely on the numerical integration of particle trajectories, requiring sufficiently resolved data in both time and space, thus, in principle, making such methods computationally expensive and time consuming [29]. Moreover, adopting a Lagrangian approach may result in a poor control over the zone of investigation, which is determined by particle motion and accumulation. For this reason, it can also be difficult to get a complete picture of the flow at specific time instants. Furthermore, the influence of particle distribution and of particle seeding schemes on quantities of interest is poorly investigated.
In contrast, Eulerian-based approaches usually simplify the data analysis workflow, as they can be directly applied to the output given by the main current techniques used for the investigation of cardiovascular flows (e.g., phase contrast MRI, CFD data). Moreover, they usually have a simpler definition, making their implementation easy and characterized by a reduced computational cost. More importantly, they can give a picture of the entire vector field. However, the inherent unsteady nature of the hemodynamic vector fields (e.g., velocity, WSS) can make the characterization of the dynamic evolution of the vector field features difficult with Eulerian-based approaches.
In conclusion, the theoretical background of the advanced methods of analysis of the WSS presented here and the recent findings related to their application to cardiovascular flows support their use to further elucidate the cause-effect relationships at the basis of the links between local hemodynamics and vascular disease. Based on the reported evidence about the physiological significance of the WSS topological skeleton in cardiovascular flows, its application in future studies, including longitudinal data, biological mechanism, and mechanobiology studies, is strongly encouraged and warranted.