LONGITUDINAL DISPLACEMENT IN VISCOELASTIC ARTERIES: A NOVEL FLUID-STRUCTURE INTERACTION COMPUTATIONAL MODEL, AND EXPERIMENTAL VALIDATION

Recent in vivo studies, utilizing ultrasound contour and speckle tracking methods, have identified significant longitudinal displacements of the intima-media complex, and viscoelastic arterial wall properties over a cardiac cycle. Existing computational models that use thin structure approximations of arterial walls have so far been limited to models that capture only radial wall displacements. The purpose of this work is to present a simple fluid-struture interaction (FSI) model and a stable, partitioned numerical scheme, which capture both longitudinal and radial displacements, as well as viscoelastic arterial wall properties. To test the computational model, longitudinal displacement of the common carotid artery and of the stenosed coronary arteries were compared with experimental data found in literature, showing excellent agreement. We found that, unlike radial displacement, longitudinal displacement in stenotic lesions is highly dependent on the stenotic geometry. We also showed that longitudinal displacement in atherosclerotic arteries is smaller than in healthy arteries, which is in line with the recent in vivo measurements that associate plaque burden with reduced total longitudinal wall displacement. This work presents a first step in understanding the role of longitudinal displacement in physiology and pathophysiology of arterial wall mechanics using computer simulations.

1. Introduction.Recent developments in ultrasound Contour and Speckle Tracking methods make it now possible to measure in vivo radial and longitudinal arterial wall displacements [44,14,48,45,15,1,13,40].These measurements for the first time reveal that longitudinal displacement of the intima-media complex in healthy Figure 1.Top: Longitudinal and radial displacement in a carotid artery measured using in vivo ultrasound speckle tracking method.The thin red line located at the intimal layer of the arterial wall shows the direction and magnitude of the displacement vector, showing equal magnitude in longitudinal and radial components of the displacement [15].Bottom: Comparison of carotid artery wall motion in a healthy (solid line) and diabetic (dashed line) subject, measured in vivo using ultrasound speckle tracking methods [15].particularly pronounced under adrenaline conditions during which the longitudinal displacement of the intima-media complex increases by 200%, and becomes twice the magnitude of radial displacement [1].It has also been shown that smaller longitudinal displacement is associated with atherosclerotic vessels [44], or with older, diabetic subjects [15] (see Figure 1; middle and right).While the relationship between radial displacement and cardiovascular disease has been extensively studied, the role of longitudinal displacement of the intima-media complex in cardiovascular disease has been completely unexplored, and the current manuscript presents an important step in this direction.
We present here a summary of our recent progress in modeling fluid-structure interaction between blood flow and arterial walls capturing both radial and longitudinal displacement, as well as viscoelastic arterial wall properties.In our work arterial walls are modeled as a linearly viscoelastic cylindrical Koiter shell capturing radial and longitudinal displacement, while blood flow is modeled using the Navier-Stokes equations for an incompressible, viscous fluid.The dynamics of the thin structure is fully coupled to the fluid dynamics of blood flow via the following two coupling conditions defined at the fluid-structure interface: (1) the kinematic coupling condition describing continuity of velocity between the fluid and structure (the no-slip condition), and (2) the dynamic coupling condition describing balance of contact forces (the viscoelastic force exerted by the structure onto the fluid is counter-balanced by the fluid stress acting onto the structure).Details of the model are presented in Section 2.
To numerically solve this FSI problem, we present a stable, loosely coupled partitioned numerical scheme which is easy to implement, and which is modular so that, e.g., different models of arterial walls, or different fluid solvers, can be easily substituted into the FSI solver by changing the appropriate fluid or structure sub-module.Moreover, additional modules, such as those modeling, e.g., plaque formation, can be easily added to the scheme.This scheme, called the kinematically coupled β-scheme, was recently introduced in [7].The scheme is based on a time-splitting approach, designed in such a way that the accuracy of the scheme is comparable to that of monolithic schemes, such as that one proposed in [5].The motion of the fluid domain is captured by the Arbitray Lagrangian-Eulerian (ALE) method.Details of the scheme are presented in Section 2.1.
In the present manuscript we use this computational model to study FSI in the common carotid artery and in a stenosed coronary artery under physiological conditions.Numerical results are compared with measurement showing excellent agreement.Additionally, new results related to the behavior of longitudinal displacement are obtained.More precisely, we show that, unlike radial displacement, longitudinal displacement in stenotic lesions is highly dependent on the stenotic geometry.In particular, we show that in Type 3 stenotic geometry presented in Section 3.2, the magnitude of longitudinal displacement is largest, which may be associated with higher incidence of plaque rupture.We also show that in Type 2 stenotic geometry exhibiting a sharp distal angle, the post-stenotic recirculation zone is largest.Based on the results in Yoganathan et al. [50], we conclude that Type 2 stenotic lesion may have a potential for higher incidence of lesion propagation.We also show that longitudinal displacement in atherosclerotic arteries is smaller than in healthy arteries, which is in line with the recent in vivo measurements that associate plaque burden with reduced total longitudinal wall displacement [45].Details of the comparison between our numerical results and experimental measurements are presented in Section 3.

2.
The Fluid-Structure Interaction Model.We will be assuming that the arterial walls are homogeneous, isotropic, and that the average thickness of the arterial walls, h, is small with respect to the length of the arterial segment L, i.e., h << L. See Figure 2.Although arterial walls consist of several layers, each with different mechanical characteristics [22,27,28], capturing the mechanics of each layer in the coupled fluid-structure interaction problem is still computationally prohibitive.Due to the computational complexity of the underlying fluid-structure interaction problem, it is common practice in hemodynamics to "average" the mechanical properties of the multi-layered arterial wall structures.We remark here, however, that the computational scheme, presented in this manuscript, is particularly suitable for modeling the multi-layered arterial wall structure due to its modularity and simplicity.This is a topic of our current research in progress.Under the above-mentioned assumptions, the mechanical properties of arterial walls can be modeled using thin structure models, such as the thin shell equations.Indeed, in the current manuscript the elasto-dynamics of the arterial walls will be modeled by the Kelvin-Voigt linearly viscoelastic Koiter shell equations, capturing both radial and longitudinal displacement.It has been shown that approximating the viscoelasticity of arterial walls by the Kelvin-Voigt model, correlates well with experimental measurements [2,10].
To derive the linearly viscoelastic Koiter shell model that includes both longitudinal and radial components of the displacement, consider a cylindrical shell of thickness h, length L, and reference radius of the middle surface equal to R (see Figure 2).This reference configuration of the thin cylindrical shell is denoted by In addition to assumption made above, we will be assuming that the load exerted onto the shell is axially symmetric, leading to the axially symmetric displacements, so that the displacement in the θ−direction is zero, and nothing in the problem depends on θ.Thus, displacement η(z, t) will have two components, the axial component η z (z, t) and the radial component η r (z, t), where these displacements refer to the displacement of the shell's middle surface.To model structural viscosity and time-dependent motion of the shell, the displacement functions depend on both space and time.The derivative with respect to the spatial variable will be denoted by η , while the derivative with respect to the temporal variable by η.
In contrast with other shell models, the Koiter shell model accounts for the contributions of both stretching and bending (flexure).The stretching of the middle surface is measured by the change of metric tensor, while flexure is measured by the change of curvature tensor, given in cylindrical coordinates, respectively, by [12] The total energy of the Koiter shell will account for the contributions of both the stretching and bending energy.We extend this model by adding the contributions due to viscous effects by employing the Kelvin-Voigt model.In this model, the total stress is linearly proportional to strain and to the time-derivative of strain.Therefore, for a linearly viscoelastic Koiter shell model the internal (stretching) force is given by and the bending moment where A is the elasticity tensor [12] defined by: and B is the viscosity tensor given by: Here, A c = 1 0 0 R 2 is the first fundamental form of the middle surface, is its inverse, and • denotes the scalar product Parameter E is the Young's modulus and σ is the Poisson's ratio, while E v and σ v correspond to the viscous counterparts of E and σ.
To define the weak formulation, which determines the energy of the viscoelastic Koiter shell, introduce the function space The weak formulation is given by: for each t > 0, find η(•, t) ∈ V c such that where ρ s denotes the volume density of the shell.Components of the forcing term f = (f z , f r ) T are the surface densities in the reference configuration of the axial and radial force.After integration by parts, the corresponding dynamic equilibrium equations in differential form can be written as follows: where and We consider a clamped Koiter shell problem with the boundary conditions given by The location of the deformed cylinder wall with respect to the reference configuration, at time t, is given by Γ(t) = {(z + η z (z, t), R + η r (z, t)) | z ∈ (0, L)} for t ∈ (0, T ).See Figure 4. Its location is not known a priori and is one of the unknowns in the problem.We will be assuming that the points (z+η z (z, t), R+η r (z, t)) on Γ(t) define a Lipschitz-continuous function (in the Eulerian framework) g(• ; t) : (0, L) → R, g(• ; t) : z → g(z; t), ∀t ∈ (0, T ), so that the fluid domain (vessel lumen) can be written as See Figure 4. Thus, the entire boundary of the fluid domain, ∂Ω(t), consists of the inlet section Γ in , the outlet section Γ out , the symmetry boundary (r = 0) Γ 0 , and the lateral boundary Γ(t), so that (see Figure 4): Although everything in this manuscript will be presented in 2 space dimensions, the main ideas directly extend to the 3-dimensional case.We are interested in simulating a pressure-driven flow through a deformable cylinder with full, two-way coupling between the fluid flow and the motion of the lateral boundary.As is often the case in such problems, without loss of generality, we will simulate only the upper half of the cylinder supplemented by the symmetry boundary condition at the axis of symmetry of the cylinder.Thus, the reference fluid domain, i.e., the reference configuration of the cylinder, is given by Ω := {(z, r)| 0 < z < L, 0 < r < R}.The lateral boundary of Ω will be denoted by The flow of blood in medium-to-large arteries is typically modeled by the Navier-Stokes equations for a viscous, incompressible fluid where u = (u z , u r ) is the fluid velocity, p is the fluid pressure, ρ f is the fluid density, and σ is the fluid stress tensor.For a Newtonian fluid the stress tensor is given by σ = −pI + 2µD(u), where µ is the fluid viscosity and D(u) = (∇u + (∇u) τ )/2 is the rate-of-strain tensor.
For the purposes of this work, at the inlet and outlet boundary the normal stress boundary conditions are prescribed: where n in /n out are the outward normals to the inlet/outlet boundaries, respectively.These boundary conditions are common in blood flow modeling [4].At the bottom boundary r = 0 the symmetry conditions are imposed: Initially, the fluid and the structure are assumed to be at rest, with zero displacement from the reference configuration u = 0, η = 0, ∂η ∂t = 0.The fluid and structure are coupled at the fluid-structure interface Γ(t) via: • The kinematic coupling condition describing continuity of velocity at the wall (the no-slip condition) • The dynamic coupling condition describing balance of contact forces at the wall where denotes the Jacobian of the transformation from Eulerian to Lagrangian framework, and σn denotes the normal fluid stress on the reference domain Ω = (0, L) × (0, R).Here e z = (1, 0) and e r = (0, 1) are the standard unit basis vectors, and n is the outward normal to the deformed domain.
2.1.The numerical scheme.To solve the fluid-structure interaction problem ( 8)-( 20) numerically, we propose to use the kinematically coupled β-scheme, recently introduced in [7].The scheme is based on a time-splitting approach known as the Lie splitting [23], described below.
2.1.1.The Lie scheme.To apply the Lie splitting scheme, the evolution problem ( 8)-( 20) must first be written as a first-order system in time: The Lie scheme consists of splitting the full problem into P sub-problems, each defined by the operator A i , i = 1, ..., P .As the original problem is discretized in time with the time step t > 0, so that t n = n t, the Lie splitting scheme consist of solving a series of problems ∂φi ∂t + A i (φ i ) = 0, for i = 1, ..., P , each defined over the entire time interval (t n , t n+1 ), but with the initial data for the i th problem given by the solution of the (i − 1) st problem at t n+1 .More precisely, set φ 0 = φ 0 .Then, for n ≥ 0 compute φ n+1 by solving and then set φ n+i/P = φ i (t n+1 ), for i = 1, . . ..I.
Roughly speaking, for our problem ( 8)-( 20), we will have the following three subproblems which will each be solved over the time interval (t n , t n+1 ), with the initial data provided by the solution of the previous sub-problem evaluated at time t n+1 : 1.The time-dependent Stokes sub-problem coupled with structure inertia on Γ via normal stress; 2. The advection sub-problem; 3. The structure elastodynamics sub-problem.This method is first-order accurate in time, and second-order accurate in space [24].To increase the accuracy in time to second-order, a symmetrization of the scheme can be performed.In [7] it was shown that the accuracy in time of the scheme presented bellow, is comparable to that of monolithic solvers, such as that one presented in [5].

2.1.2.
The first-order system in ALE framework.To deal with the motion of the fluid domain we adopt an Arbitrary Lagrangian-Eulerian (ALE) approach [26,19,33].In the context of finite element method approximation of moving-boundary problems, ALE method deals efficiently with the deformation of the mesh, especially near the interface between the fluid and the structure, and with the issues related to the approximation of time-derivatives on fluid domains that depend on time.
An ALE approach is based on introducing a family of (arbitrary, invertible, smooth) mappings A t defined on a single, fixed, reference domain Ω such that, for each t ∈ (t 0 , T ), A t maps the reference domain Ω = (0, L) × (0, R) into the current domain Ω(t).In our approach, we define A t to be the harmonic extension of the mapping that maps the boundary of Ω to the boundary of Ω(t) for a given time t.
More precisely, we define A t to be the solution of the Laplace's equation: Differentiation with respect to time, after using the chain rule, gives ∂f ∂t x = ∂f ∂t + w • ∇f, where w denotes domain velocity given by Notation ∂f ∂t x means that the time-derivative is calculated at the reference domain Ω.
Using this knowledge, we can now write system ( 8)-( 20) in the ALE framework.At the same time, we will write this system in first-order form so that we can immediately apply the Lie splitting.To write system ( 8)-( 20) in first-order form we will utilize the kinematic coupling condition (18) and express the second-order time derivative of structure displacement in terms of the first-order derivative of the trace of the fluid velocity at the structure boundary û| Γ = û(ẑ, R, t).
Written in the ALE framework, our problem now reads: with the kinematic coupling condition on Γ: and the dynamic coupling condition on Γ (incorporating condition (29)): where J is the Jacobian given in (21).The following boundary conditions on Γ in ∪ Γ out ∪ Γ 0 are enforced: At time t = 0 the following initial conditions are prescribed: , where β is a number between 0 and 1, 0 ≤ β ≤ 1, with β = 0 corresponding to the splitting introduced by Guidoboni et at.[24].Part I of the fluid stress will be used to load the viscous part of the structure, treated together with the fluid in Step 1 of the splitting (the time-dependent Stokes problem).Part II of the fluid stress (the pressure) will be used to load the pure elastodynamics problem for the structure in Step 3. Thus, parameter β is used to distribute the fluid pressure between the fluid and structure sub-problems.Numerical investigation in [7] showed that the change in β affects the accuracy of the scheme, but not stability.In fact, it was proved in [11] that our scheme is unconditionally stable for all 0 ≤ β ≤ 1.For the problem studied in [7], the value of β = 1 provided the highest accuracy.For this reason, β = 1 will be used in the splitting presented in this manuscript.More precisely, the splitting scheme is given by the following.
Step 1.The time-dependent Stokes problem is solved on a fixed domain Ω(t n ) with the boundary condition on Γ(t n ) which couples the normal fluid stress with structure inertia, and with the viscous energy of the structure (when it is non-zero).This can be done by using the kinematic coupling (implicitly) to express the timederivatives of structure displacement in terms of the trace of fluid velocity on Γ(t n ).The problem reads as follows: Find u, p and η, with û(x, t) = u(A t (x), t) such that for t ∈ (t n , t n+1 ), with p n and η n obtained at the previous time step, the following is satisfied: with the following lateral boundary conditions on Γ(t n ): where The boundary conditions on Γ in ∪ Γ out ∪ Γ 0 are given by: and initial conditions by Then set u n+1/3 = u(t n+1 ), η n+1/3 = η(t n+1 ), p n+1 = p(t n+1 ).
Notice how in Part I of the fluid stress, β pn is taken from the previous time step, while the fluid stress σ is treated implicitly, at the current time step.The Jacobian J, and the normal n to the lateral boundary, are taken from the previous time step.
Step 3: This step involves solving an elastodynamics problem for the location of the deformable boundary by involving the remaining, elastic part of the structure, which is loaded by Part II (pressure) of the normal fluid stress.The pressure is obtained from the Stokes problem (Step 1), at the "current" time step.Additionally, the fluid and structure communicate via the kinematic coupling condition which is used as the initial velocity for this structural problem.The problem reads: Find û and η, with pressure p n+1 computed in Step 1, and Jacobian and normal n n obtained from the previous time step, such that for t ∈ (t n , t n+1 ) with boundary conditions: and initial conditions: Then set u n+1 = u(t n+1 ), η n+1 = η(t n+1 ).Do t n = t n+1 and return to Step 1.
Remark 1.We remark here that no sub-iterations between the fluid and structure sub-solvers are needed for stability.In fact, it was shown [11] that this scheme is unconditionally stable even when the density of structure is roughly equal to the fluid density, which is the case in hemodynamics applications.
Remark 2. We emphasize, again, that the main feature of this stable, loosely coupled scheme, is the approach presented in Step 1.In Step 1 a portion of the "coupled" FSI sub-problem is solved on a fixed domain, where the coupling is done through the structure terms involving structure velocity, i.e., the inertia term, and viscous structure effects.The structure velocity is expressed, via the kinematic coupling condition, in terms of fluid velocity (trace) on Γ(t n ), which enables writing the entire problem as a fluid problem on Ω(t n ), where the normal fluid stress on Γ(t n ) is balanced (in an implicit way) by a portion of the structure load.This way the kinematic coupling condition and a (partial) dynamic coupling condition are implicitly enforced in the fluid solver already at the first step of the scheme.This tight coupling between the fluid and structure inertia, and therefore, between the fluid and structure kinetic energy, is responsible for the stability of this scheme.
A block diagram summarizing the main steps of the scheme is shown in Figure 5. 3. Results.In this section we show that our computational model gives rise to the physiologically reasonable solutions by studying two examples of problems for which there exist data obtained from the measurements of both radial and longitudinal displacement.The first example concerns a healthy common carotid artery (CCA), while the second example concerns an atherosclerotic coronary artery.New results related to the influence of the geometry of stenotic lesion on the magnitude of longitudinal displacement will be shown.
We mention here that our computational model was tested against other numerical solvers on a benchmark problem for FSI in hemodynamics, for which only radial displacement was compared with existing simulations, showing excellent agreement and performance (see Bukac et al. [7]).Our method was shown to be first-order accurate in time, and second-order accurate in space (see Guidoboni et al. [24]).The accuracy of the method was shown to be comparable to that of the monolithic scheme by Badia et al. [5,7].
In the current manuscript, for the first time we compare solutions of the presented computational model with experimental measurements, which we present next.
3.1.The common carotid arteries (CCA).Parameter Values.Left and right common carotid arteries in human subjects differ significantly in length and in their mode of origin.The study in Ribeiro et al. [41] measured average length of the left and right CCA in 46 male cadavers.The measured length of the right CCA was 9.6 ± 0.1 cm while the left CCA measured 12.1 ± 0.2 cm.The diameter of CCA slightly differs in males and females, and ranges between 6.5 ± 0.99 cm in males and 5.97 ± 0.9 cm in females [31,30].Measurements in Bussy et al. [9] reported wall thickness of 0.0582 ± 0.0139 cm.
Our choice of parameters lies in the bounds given above, and is summarized in Table 1.The viscoelastic constants that appear in our structure equations are equal 1.055 Fluid viscosity µ (g/(cm s)) 0.04 Wall density ρ s (g/cm 3 ) 1.055 Wall thickness h (cm) 0.07 Young's mod.E(dynes/cm 2 ) 2 × 10 6 Poisson's ratio σ 0.5 Table 1.Geometry, fluid and structure parameters for common carotid artery.

MARTINA BUKA Č AND SUN ČICA ČANI Ć
to This choice of viscoelastic parameters is within the range of measured viscous moduli of blood vessels reported in Armentano et al. [3].
Pressure Data.We study blood flow driven by the inlet and outlet pressure data shown in Figure 6.The shape of the pressure wave is the same as in Warriner et al. [48] while the pressure drop is scaled by the factor 0.9 to recover physiologically reasonable blood velocity and Reynolds number.Namely, several experimental studies have shown that longitudinal velocity in healthy common carotid artery is usually less than 100 cm/s [42,6,17,32] giving rise to the local Reynolds number which is less than 1500.This is smaller than the results in Warriner et al. [48] which recover blood flow conditions corresponding to the local Reynold's number over 2400.Indeed, as we shall see below, our computed velocity is in very good agreement with the Doppler velocity data reported by Weinberg et al. [49].
Fluid Velocity.We compared the velocity computed using the data listed above, with Doppler velocity measurements of the left common carotid artery (CCA), reported by Weinberg et al. [49].Figure 7, bottom right, shows the magnitude of velocity at the center of the vessel over a cardiac cycle calculated using our computational model.This should be compared with the Doppler velocity graph, shown in Figure 7, top right.One can see that the graphs agree both quantitatively and qualitatively.The slight difference in the morphology of the graphs is due to the fact that we did not have the exact pressure drop data for this particular velocity measurement reported in [49].However, we see that the maximum velocity obtained using our computational model is around 95 cm/s, and the minimum around 30 cm/s.This is in agreement with several results showing experimentally measured CCA velocities [42,6,17,32].In particular, the measurement shown in Figure 7 shows the maximum velocity of 101.1 cm/s, and the minimum velocity of 28.4 cm/s.This shows the difference of only 6% for the maximum, and 3% for the minimum velocity when compared with our computational results.placement decreases with age, and usually varies between 0.1 mm and 0.38 mm for the CCAs, i.e. between 3% and 13% of the vessel radius.Recent experimental studies obtained using B-mode ultrasound speckle tracking method and/or B-mode ultrasound velocity vector imaging, report longitudinal displacement of the same magnitude as radial displacement [44,13,40,45].In Table 2 we report the measurements of the diameter change and total longitudinal displacement of CCA, as reported in Svedlund et al. [44].These are compared to the results obtained using our computational model, showing that our results fall precisely within the confidence interval of the measurements.Notice also that the computed maximal radial displacement is around 6%, which is well within the normal range [9,17].
Figure 8 shows the corresponding numerically computed radial and longitudinal displacement curves at the midpoint of the vessel wall, over one cardiac cycle.
Viscoelasticity.The viscoelastic effects are visible in the stress-strain relationship of the arterial wall, which exhibits hysteresis.In our simulations we used the viscosity parameters for the vessel wall, listed in (38), which are within the range of parameters reported in Armentano et al. [3].We calculated the hysteresis behavior in terms of diameter vs. pressure, and obtain the curve depicted in Figure 9 (left).
To quantify the hysteresis behavior one can calculate the Energy Dissipation Ratio (EDR), which is a measure of the area inside the diameter-pressure loop relative to the measure of areas inside and under the loop.More precisely, if A 1 is the area inside the hysteresis loop, and A 2 is the area below the hysteresis loop, then EDR is defined to be EDR = A 1 /(A 1 + A 2 ) × 100%.Walls with higher viscoelasticity have larger area inside the loop, resulting in higher EDR.In our simulations EDR is 8.5%.This is comparable to the results in Warriner et al. [48] which show EDR of around 7.8% for younger subjects.

Coronary artery stenosis.
A stenosed coronary artery with mild stenosis (60%) is considered.Three different geometries of stenotic lesions (see Figure 10) were studied for radial and longitudinal displacement: symmetric geometry (type 1), a geometry with a small inflow angle and a sharp outflow angle (type 2), and a geometry with a sharp inflow angle and a smooth outflow angle (type 3).Parameter Values.Experimental study in Mokhtari et al. [38] reported Young's modulus of a healthy coronary artery of 3.84 ± 0.389 × 10 6 dynes/cm 2 , while in mild stenosis it increased up to 5.02 ± 0.07 × 10 6 dynes/cm 2 .In our numerical example, the Young's modulus is a function of position that changes between those two values  3. Geometry, fluid and structure parameters for coronary artery.
depending on the geometry.The reference radius for a healthy coronary artery is 0.18 cm, as measured in Dodge et al. [18].The reference radius in stenotic lesions was calculated based on the following formulas: Since the reference radius and Young's modulus now depends on the position, we derived an appropriate model for the structure from the Koiter shell model.Details on the derivation of the model (with only radial displacement) can be found in Tambaca et al. [46].The values of the parameters used in the simulations are given in Table 3.
Pressure data.The inlet and outlet pressure were taken from the measurements in Marques et al. [36] where a trans-stenotic pressure gradient in coronary arteries was recorded.The data is shown in Figure 11.We utilized those values of the pressure gradient as the inlet and outlet data in our simulations.Velocity.We calculated the fluid velocity wave form and streamlines for the 3 types of stenotic lesions.The numerically calculated velocity wave form over one cardiac cycle is shown in Figure 12, bottom.The numerically computed peak systolic velocity is around 25 cm/s in the area proximal to stenosis, and 60 cm/s in the stenotic region.This is in very good agreement with the measurements reported in Hozumi et al. [25], in Johnson et al. [29], and in Marques et al. [36].The results reported in Marques et al. [36] are shown in Figure 12, top.They report the measured peak systolic velocity in the area proximal to stenosis around 28 cm/s (compare with 25 cm/s, i.e., 10% difference with numerical simulation), and in the stenotic lesion around 70 cm/s (compare with 60cm/s i.e., 14% difference with numerical simulation).
Figure 12 also shows a very good agreement between the morphology of the measured and numerically calculated velocity wave curves over one cardiac cycle.We emphasize here that, unfortunately, the exact inlet/outlet pressure was not given in Marques et al. [36], only the pressure gradient over a stenotic lesion.The inlet and outlet pressure data that we used in our simulations is shown in Figure 11.The corresponding trans-stenotic pressure gradient for this data is the same as the pressure gradient reported in Marques et al. [36].
Notice that our results indicate that the velocity wave form over one cardiac cycle does not change significantly between the three types of stenotic lesions, however, it significantly increases in the stenotic region, as expected.
Figure 13 shows the numerically calculated streamlines for the 3 different stenotic geometries, recorded near the peak pressure gradient, at t = 0.35 s.One can see that for the type 2 stenotic geometry with a sharp distal (outflow) stenotic angle, the post-stenotic region exhibits a larger and more pronounced recirculation zone and stagnation points, which can be a pre-cursor for a propagation of stenotic lesion.This is in line with the observations presented in [20,50], where it was remarked that in recirculation zones the shear stress is high, and the platelets trapped in the zone experience a longer exposure time, which have been associated with platelet aggregation and activation.
Radial and Longitudinal Displacement.We simulated radial and longitudinal displacement of the arterial wall proximal to stenosis and in the stenosed area for the three different types of stenotic lesions, see Figure 14.By comparing the radial displacement proximal to stenosis (Figure 14 top left) and in the stenosed area (Figure 14 bottom left), one can immediately observe that the radial displacement in the stenotic region is two orders of magnitude smaller than in the region proximal to stenosis (5×10 −5 vs. 10 −3 ).This is due to the local wall stiffening which is associated with atherosclerosis, and is captured in our model by the increase in the Young's modulus at the location of stenosis.Another interesting observation is that radial displacement does not show significant differences between the three types of stenotic geometries.The most interesting new observation, however, is  that the longitudinal displacement is very different for the three different stenotic geometries.In particular, the plots on the right in Figure 14, show that longitudinal displacement is largest in a coronary artery with the type 3 stenotic lesion.This may indicate that type 3 geometry may be associated with higher incidence of plaque rupture.Further research in this area is needed to make that correlation.
Finally, we compared the magnitude of longitudinal displacement in a healthy coronary artery with a coronary artery suffering from atherosclerosis (type 1 stenotic lesion).See Figure 15.The goal was to see if our simulations can predict the result presented in the recent work by Svedlund at al. [45] in which it was demonstrated, using velocity vector imaging (Vevostrain), that plaque burden in atherosclerotic arteries is associated with lower total longitudinal wall motion.It was hypothesized that the main reason for this observation is the reduced longitudinal tensile stress which is typical in local wall thickening due to the higher hemodynamics stress.This phenomenon should be captured in our model by the increase in the Young's modulus which affects both the magnitude of radial as well as longitudinal wall motion, which can be seen by the form of the coefficients in the Koiter shell model ( 8), (9), in which both the radial and longitudinal stress coefficients depend on the Young's (stiffness) modulus of the structure.Indeed, Figure 15 shows exactly that: the maximum longitudinal displacement in a healthy artery is significantly bigger (twice the size) than the longitudinal displacement in an atherosclerotic artery.Viscous Dissipation.We conclude this manuscript by showing the viscous dissipation for the three different stenotic geometries.Viscous dissipation measures, among other things, flow disturbances caused by stenoses.Larger flow disturbances lead to higher values of the dissipated energy [35,34,16].Figure 16 shows the values of viscous dissipation, integrated in time over one cardiac cycle.The viscous dissipation 2µ|D(u)| 2 was integrated over the cross-section of the vessel S i , for each point z i = i h on the horizontal axis, where i = 0, ..., 340, and h = 0.0176 cm.Thus, the three graphs in Figure 16 show the values of s(z) = S(z) 2µ|D(u)| 2 dy, for each of the three stenotic geometries.It is interesting to notice that the stenotic geometry of type 1, for which the proximal and distal angles of the stenotic lesion are smaller than the distal angle in type 2 geometry, and smaller than the proximal angle in type 3 geometry, exhibits smallest viscous energy dissipation.In fact, we have calculated the total viscous energy dissipation for each of the 3 cases by integrating the values of viscous energy dissipation over space and time, to obtain that the total viscous energy dissipation equals 1322.22 g/(cm • s 2 ) for type 1 geometry, 1817.21g/(cm • s 2 ) for type 2 geometry, and 1921.04 g/(cm • s 2 ) for type 3 geometry.Thus, the total viscous energy dissipation for type 1 geometry is 32% smaller than the total viscous energy dissipation for type 3 geometry.Therefore, even though degree of stenosis in all three cases is the same, the stenosis with the smallest proximal and distal angles exhibits the smallest viscous energy dissipation.Notice that even though the vortices in type 3 geometry seem to be the smallest, the total loss of viscous energy is still large in this case, due to the large losses associated with the flow impinging the entrance to the stenotic lesion with a relatively large incident angle.

4.
Discussion.This work presents a first step in the understanding of the role of longitudinal displacement in the physiology and pathophysiology of the human cardiovascular system using computational modeling.The dynamics of longitudinal (and radial) displacement is captured by the linearly viscoelastic Koiter shell model, which is fully coupled to the Navier-Stokes equations for an incompressible, viscous fluid modeling blood flow.A novel partitioned, loosely coupled scheme was presented for the simulation of the coupled fluid-structure interaction problem.
Physiologically relevant examples were considered for which the results of computer simulations were compared with experimental data showing excellent agreement.New results related to the behavior of longitudinal displacement were obtained.More precisely, it was shown that, unlike radial displacement, longitudinal displacement in stenotic lesions is highly dependent on the stenotic geometry.In particular, it was shown that in type 3 stenotic geometry, the magnitude of longitudinal displacement is largest, which may be associated with higher incidence of plaque rupture.It was also shown that longitudinal displacement in atherosclerotic arteries is smaller than in healthy arteries, which is in line with the recent in vivo measurements that associate plaque burden with reduced total longitudinal wall displacement [45].

Figure 2 .
Figure 2. A sketch of the arterial wall.

Figure 3 .
Figure 3. Left: Cylindrical shell in reference configuration with middle surface radius R and shell thickness h.Right: Deformed shell.

) 2 . 1 . 3 .
Details of the operator-splitting scheme.To administer the modified Lie splitting mentioned of Section 2.1.1,we first split the fluid stress σn into two parts, Part I and Part II, as follows: σn = σn + β pn

Figure 5 .
Figure 5.A block diagram showing the main steps of the numerical scheme.

Figure 8 .
Figure 8. Numerically calculated radial and longitudinal displacements over one cardiac cycle.

Figure 10 .
Figure 10.A sketch of three different stenotic geometries used in the study.

Figure 13 .
Figure 13.Comparison of the velocity and streamlines for the three different stenotic geometries at t = 0.35 s.

Figure 14 .
Figure 14.Radial and longitudinal displacement of the vessel wall proximal to stenosis (top) and in the stenosed area (bottom).

Figure 15 .
Figure 15.Longitudinal displacement in a healthy and stenosed coronary artery.

Figure 16 .
Figure 16.Viscous dissipation for the 3 types of stenotic geometries, shown in Figure13.