Geometry and Flow Properties Affect the Phase Shift between Pressure and Shear Stress Waves in Blood Vessels

: The phase shift between pressure and wall shear stress (WSS) has been associated with vascular diseases such as atherosclerosis and aneurysms. The present study aims to understand the effects of geometry and ﬂow properties on the phase shift under the stiff wall assumption, using an immersed-boundary-lattice-Boltzmann method. For pulsatile ﬂow in a straight pipe, the phase shift is known to increase with the Womersley number, but is independent of the ﬂow speed (or the Reynolds number). For a complex geometry, such as a curved pipe, however, we ﬁnd that the phase shift develops a strong dependence on the geometry and Reynolds number. We observed that the phase shift at the inner bend of the curved vessel and in the aneurysm dome is larger than that in a straight pipe. Moreover, the geometry affects the connection between the phase shift and other WSS-related metrics, such as time-averaged WSS (TAWSS). For straight and curved blood vessels, the phase shift behaves qualitatively similarly to and can thus be represented by the TAWSS, which is a widely used hemodynamic index. However, these observables signiﬁcantly differ in other geometries, such as in aneurysms. In such cases, one needs to consider the phase shift as an independent quantity that may carry additional valuable information compared to well-established metrics.

Hemodynamic forces, such as pressure and wall shear stress (WSS), play an important role in vascular physiology and pathology [6][7][8]. The interaction between tensile and shear stresses has been widely reported to play a crucial role in endothelial behavior [9,10]. The flow-induced WSS is time varying due to the pulsatile nature of blood flow. The phase shift (temporal delay) between the pressure and WSS on the endothelium is a primary parameter that characterizes the interaction between these two forces [11] and is closely associated with vascular diseases [12,13].
It has been suggested that the phase shift is important in vascular regulation and remodeling [14] and in predicting sites prone to atherosclerosis in the coronary artery [11], the carotid artery [15], and arterial stenosis [16]. It has been proposed by several independent researchers that a large phase shift occurs in regions where atherosclerotic lesions are prevalent [16][17][18][19][20]. Moreover, Maul et al. [21] reported that the magnitude and/or frequency of these mechanical stimuli affect the morphology, growth, and differentiation of mesenchymal stem cells in the vasculature. Owatverot et al. [22] performed an analysis of concurrent tensile and shear stresses and highlighted their influence on the behaviors of endothelial and stem cells. The phase shift has also been shown to significantly affect the strain energy density of the endothelial cell where mechanotransduction is prevalent [23]. Shojaei et al. [24] reported that the phase shift regulates the differentiation of human-adipose-derived stem cells toward the endothelial phenotype.
It has been further noted that the phase shift influences proatherogenic responses in endothelial cells and is an important parameter characterizing arterial susceptibility to disease [25,26]. These observations are consistent with the understanding at the cellular or molecular level: the endothelial monolayer senses the phasic change between pulsating pressure and WSS and alters its cellular actions, in particular nitric oxide production [19,27]. Atherosclerotic plaques have been reported to be also closely associated with a low timeaveraged WSS (TAWSS) [28][29][30] and a high temporal WSS gradient (TWSSG) [31][32][33].
As to the underlying physical mechanism for the phase shift, it is noteworthy that pressure fluctuations propagate with the speed of sound. This process is fast compared to variations of the fluid velocity field and of its spatial gradients, the latter determining the stress tensor. As a consequence, one expects a time delay and thus a phase shift between pressure and shear stress signals. While this issue can be addressed by analytic means for a pulsatile flow in a straight tube (see the Womersley solution below), numerical modeling is necessary to address this issue in more complex geometries.
Indeed, little is known about the effects of complex vessel geometry and fluid flow properties on the phase shift and on the possible correlation among the WSS-based quantities, such as the phase shift, the TAWSS, and the TWSSG. There is probably a coupled causality: geometry alters flow (short-term effect), and flow changes geometry (long-term effect). In the present study, we aim to investigate the first case and assess the hypothesis that geometry alone can play a crucial role in the phase shift. In principle, the pulsatile deformation of an artery (or aneurysm) may also have some effects on the time delay. As a first step, however, the focus of this work is to understand the roles the geometric and flow parameters play under the rigid (or stiff) [34,35] wall assumption.
In this work, we employed lattice-Boltzmann simulations of flow in differently shaped blood vessels ( Figure 1) to study the change of the phase shift with geometry, the Womersley number (flow pulsatility), and the Reynolds number (flow speed). We found that geometry has a strong effect on the phase shift. While the phase shift is closely correlated with other observables, such as the TAWSS, in straight or slightly curved blood vessels, it was observed to become unpredictable within an aneurysm, which might accelerate the development of vulnerable regions and eventually the rupture of the blood vessel.
The paper is organized as follows. Relevant physical aspects are introduced in Section 2. The numerical methods and benchmark tests are presented in Sections 3 and 4. The simulation results of sinusoidal flow in differently shaped domains are presented and discussed in Section 5, with a particular focus on the geometrical effect. Section 6 provides a summary of the work with an outlook for future studies. (b,c) represent a curved blood vessel without and with a side-wall aneurysm, respectively. The membrane that defines the blood vessel has zero thickness and is represented by a set of vertices and flat triangular facets.

Physical Model
Throughout this study, it was assumed that the fluid (blood) can be represented as an incompressible fluid with uniform density. The continuum assumption is justified by the large diameter of the blood vessel considered in this study (a few millimeters up to centimeters) as compared to the typical size of the particulate objects (red blood cells, platelets, and white blood cells, which typically lie in the range 2-20 µm). Thus, the particulate properties [36,37] of the blood were not considered here. We thus sought to solve the underlying Navier-Stokes equations in the incompressible limit. It is noteworthy that the methodology implemented here-the lattice Boltzmann method (LBM)-allows for a small, but finite numerical compressibility. It turns out that this feature is not necessarily a drawback. Rather, if used adequately-e.g., by introducing multiple relaxation times (MRTs) within LBM; see Section 3-it may lead to better stability of the numerical scheme.

Physical Ingredients
In principle, the circumferential stretch of a blood vessel is dominated by the socalled transmural pressure, which represents the pressure difference on both sides of the vessel wall. The external (or intracranial) pressure is normally one order of magnitude smaller than the intra-artery pressure and is often neglected in most existing models [38][39][40][41]. Hence, the temporal variation of the transmural pressure (or circumferential stretch) mainly depends on the pressure inside the artery or aneurysm.
The pressure pulse generated by the heart is a composite wave consisting of multiple harmonic components. Each harmonic component is a sine or cosine wave with different amplitude and frequency. In this work, we used a single sine curve to represent the input signals in order to focus on fundamental physical aspects.
Without loss of generality, this work: (1) used a toy model and focused on a short segment (patch) of the blood circulatory system; (2) assumed that the waveform of pressure at the inlet and outlet is known and sinusoidal; any change in the inlet/outlet pressure was treated as an input to the simulation; (3) focused on WSS since the pressure signal has a long wavelength and its traveling time is negligibly short (see Appendix A for more details).
The phase shift between the pressure, P, and WSS, σ w , is characterized by the phase shift, ∆ϕ, in degrees (as adopted in [42,43]): where t(P max ) and t(σ max w ), respectively, indicate the peak pressure and WSS times within one cardiac cycle of duration T (Figure 2). The phase shift is positive when the pressure signal reaches its maximum before the WSS signal and negative if the pressure lags behind the WSS.
Here, t(σ max w )/T (blue bar) and t(P max )/T (red bar) represent the times of maximum WSS and pressure within one period T, respectively. This kind of behavior is examined at each surface point in the vessel, as each point may have a different time evolution of WSS.
As discussed in Appendix A, it is reasonable to assume that the maximum pressure time at any point of the short vessel considered in this work is equal to the peak time of the input pressure signal, i.e., t(P max ) = T/4 for the entire vessel since the pressure signal in our model is a sine wave.
Besides the phase shift, the temporal distribution of WSS is also quantified by the timeaveraged wall shear stress (TAWSS) and the temporal WSS gradient (TWSSG). The TAWSS is determined by integrating the WSS magnitude at each point x over one period [44]: The TWSSG quantifies the average rate of change in WSS over half a cardiac cycle [32,33], where σ max w (x) and σ min w (x) indicate the maximum and minimum values of WSS at x and within the considered time interval.

Dimensionless Groups
The dimensionless numbers that are relevant to this work can be defined as ratios of the physical time scales mentioned in Appendix A: the Reynolds number, which describes the ratio of inertial to viscous interactions, Re = UD/ν = (2L/R)t ν /t u , where R and L, respectively, indicate the radius and length of the vessel (L/R = 12 and diameter D = 2R), ν is the kinematic viscosity, U indicates the mean flow velocity, and t ν and t u indicate the momentum diffusion and advection times, respectively; the Womersley number, which characterizes the pulsatile flow frequency in relation to viscous effects, where ω is the angular frequency (ω = 2π/T with T indicating the period). The last ratio is the Mach number, Ma ∼ t s /t u , with t s representing the sound wave (or acoustic) time, which is not of interest in the context of this work as it is always small in blood flow.

Numerical Model
In this work, the standard incompressible Navier-Stokes equations of fluid dynamics in three dimensions were solved using a standard lattice Boltzmann method on a D3Q19 lattice with the multiple-relaxation-time collision operator [45]. The coupling between the fluid and vessel wall was realized via the immersed boundary method [46]. The use of the MRT scheme within LBM allows tuning the relaxation parameters such that both the accuracy and stability of the numerical scheme are improved [45]. This feature is particularly important at high Reynolds numbers, where small fluctuations may quickly grow due to enhanced nonlinearity. The method and the underlying in-house software were tested versus semi-analytic results by Krüger et al [47]. Further applications of the software include suspensions of deformable red blood cells [36,48,49] and compliant blood vessels [50].
The effect of the mesh size on the numerical accuracy of the present model was already addressed in [47]. Based on the experience gained from this work and a more recent study [50], we found that, within given constraints on the computation time, a mean effective spatial resolution (mesh size) of approximately 0.087 mm provides very good results. As will be shown below (see Section 4), with this choice, the simulations reproduced the predicted phase shift within an error of less than 2%. The fluid domain consisted of approximately 3.7 million lattice nodes. The blood vessel (immersed structure) consisted of approximately 82,000 nodes. The discrete time step was set to dt = 6.7 × 10 −6 s; for a cycle with a duration of 0.8 s, the temporal resolution was approximately 120,000 time steps per cycle. Typically, each cycle ran for approximately 40 h in serial on a workstation with Intel® Xeon(R) CPU E5-1620 v2 @ 3.70GHz.
The limit of a quasi-rigid tube wall was achieved by assigning stiff springs on the entire wall, together with making all elastic moduli as large as possible; this is the same treatment as in [50].
In lattice-Boltzmann simulations, the numerical speed of sound is generally smaller than that in reality; it was 7.6 m/s in our simulations. Nevertheless, all physical time scales (Appendix A) that were resolved were sufficiently long compared to the numerical sound propagation time. Using the parameters of the present study, it is easy to verify that t ν = 356t s , t ω = 270t s , t u = 44t s , and t ρu = 45t s , where t ω and t ρu indicate the oscillation and inertial times, respectively. Neglecting sound propagation is thus a good approximation in the present study. Throughout this work, the peak WSS time, t(σ max w ), was extracted by curve fitting to a sinusoidal function, using simulated data within a short time interval of T/4 around each peak time. As justified in Appendix A, we set for the entire tube t(P max ) = T/4.
We used the pressure-corrected periodic boundary conditions [51] to prescribe an intended pressure drop between the two ends of the blood vessel to drive flow, together with applying periodicity to fluid velocity. The no-slip condition at the wall was recovered by the employed immersed boundary method.

Validation
The accuracy of the present model in the quantitative analysis of the phase shift (or peak shear stress time, t(σ max w )) is demonstrated as follows. For flow through a straight pipe, driven by a sinusoidal pressure gradient: the phase shift ∆ϕ between waves of WSS and pressure can be expressed as a function of the Womersley number α [52], where: with ber (·) and bei (·) representing the derivatives of Kelvin functions ber(·) and bei(·), respectively. Simulation parameters relevant for this validation are given in Table 1. We used a Python script to calculate the analytical phase shift by evaluating Equation (5). Throughout this work, we used Equation (1) to determine the phase shift from the simulations. As illustrated in Figure 3, the maximum relative error of our simulation results from the corresponding analytical solutions was less than 2% when averaging over the entire tube wall. Without such averaging, the maximum local error was 5.7%. Table 1. Simulation parameters for the benchmark tests in simulation units. The simulation box is (N x , N y , N z ) = (256, 57, 57). The fluid density is ρ = 1, and the tube radius is R = 22.9 (tube length L = 12R). We apply stiff springs on all Lagrangian nodes to mimic nearly rigid walls. The elastic shear modulus κ s , area dilation modulus κ α , and spring modulus κ sp have the same magnitude of 0.3. The bending modulus is κ b = 0.03. All other relevant parameters are listed below. P (t): pulsatile pressure gradient (note Equation (4) with P 0 indicating the amplitude); T: period; η: dynamic viscosity. For blood flow in major brain arteries, the Womersley number, α, is around 3. In validation tests, the Womersley number varies from 1.02 to 2.88 via changing the period and/or viscosity. The Reynolds number is either 5 or 179.  The analytical data (black line) are obtained via Equation (5) using a Python script. In simulations using the geometry shown in Figure 1a, we alter the viscosity ν and/or period T to control α. The simulated phase shift is averaged over the entire tube wall (red circles). We observe an excellent agreement, with the maximum numerical error being less than 2%. The inset shows the analytic result for large α, suggesting that ∆ϕ → 45 • for α → ∞.

Results and Discussion
Unless otherwise stated, the simulation parameters given in Table 2 were used; the sinusoidal pressure gradient shown in Equation (4) was used for all simulations. The corresponding Reynolds and Womersley numbers were Re = 179 and α = 2.88, respectively, which are typical for cerebral blood flow. Based on sinusoidal flow in a straight cylindrical tube, we can define three characteristic quantities: spatiotemporal average of flow velocity, U = P R 2 /8η = 0.17 m/s; time-averaged wall shear stress, TAWSS = 3P 0 R/4 = 2 Pa; temporal WSS gradient, TWSSG = P 0 R/T = 3.4 Pa/s.

Study Cases
To disentangle the effects of wall geometry and fluid flow on the phase shift between pressure and WSS, we considered two different sets of simulations:

1.
First, we focused on the effect of flow. We performed the simulations of sinusoidal flow through a curved tube without an aneurysm (Figure 1b). We varied the flow properties in two ways: In the second set of simulations, we used three different geometries ( Figure 1) and kept the flow properties identical: the average flow velocity, Reynolds, and Womersley numbers were 0.17 m/s, 179, and 2.88, respectively (see Figure 6).
Each simulation was performed for a total duration of three cardiac cycles, which is sufficient to reach a time-periodic flow field.

Varying Flow Properties in a Curved Tube without Aneurysm
The phase shift between the pressure and WSS waves, ∆ϕ, is affected by several factors. We started with the role of fluid flow properties, in particular the Womersley and Reynolds numbers, α and Re.
In a straight tube, ∆ϕ is the same on every point of the tube surface, and it increases with α, ranging from 0 • to 45 • (Figure 3). We observed that, in a curved tube, ∆ϕ becomes inhomogeneous ( Figure 4). While the average value of ∆ϕ increases with α, the variation of ∆ϕ also increases with α, and ∆ϕ can reach local values larger than 45 • for sufficiently high α. In agreement with the calculations in [53], we observed that the TAWSS and TWSSG increase with α. Additionally, our simulations showed that both TAWSS and TWSSG become increasingly heterogeneous with α. A key feature emerging with increasing α is the formation of regions opposite each other where ∆ϕ, the TAWSS, and the TWSSG assume low and large values, respectively. Following the flow entering the tube from the left in Figure 4, the flow is first forced upward at the bend, resulting in an increase of ∆ϕ, the TAWSS, and the TWSSG at the bottom of the tube and a decrease at the top. In the second half of the bend, when the flow is forced downward again, the trend reverses, and the observables assume larger value at the top of the tube and lower values at the bottom.
As shown by Equation (5), ∆ϕ in a straight tube is a function of α only. In particular, ∆ϕ does not depend on Reynolds numbers. Our simulations confirmed the independence of ∆ϕ on Re when the tube is not curved (data not shown). For a curved pipe, however, we found that Re also largely controls ∆ϕ; a higher Re leads to a larger ∆ϕ ( Figure 5). The spatial distribution of ∆ϕ is homogeneous at the lowest Re investigated (Figure 5a). It then becomes increasingly heterogeneous as Re is raised from 5 to 179 (Figure 5b,c). This spatial heterogeneity seems to saturate and even decrease upon the further increase of Re (Figure 5d). Surprisingly, ∆ϕ can even reach negative values in certain places for sufficiently large Re. Similar to Figure 4, increasing Re also leads to the formation of two pairs of oppositely located patches with locally small and large values of ∆ϕ, the TAWSS, and the TWSSG.
It is noteworthy that the variations of the simulation parameters have no impact on ∆ϕ, as long as Re and α are kept constant. For example, compared with the case shown in Figure 4a, a simulation with 3.4-times higher flow speed, viscosity, and frequency (but fixed α = 1.56 and Re = 50) gave the same distribution and values of ∆ϕ (data not shown). This result simply confirms the law of similarity, which states that dimensional parameters, such as the flow rate and period, influence the flow field (except for its absolute magnitude), only through dimensionless numbers, such as Re and α.

Curved Tube with Aneurysm
We saw that, besides the flow properties, the geometry also affects the phase shift ∆ϕ between the pressure and shear stress signals. In the case of pulsatile flow in a straight tube, ∆ϕ is the same everywhere in the tube (Figure 6a) if we neglect small numerical errors due to discretization of the geometry. However, we found that ∆ϕ varies in space for other geometries, such as the curved tube ( Figure 6b) and the curved tube with aneurysm ( Figure 6c). Figures 4 and 5 show that the effect of geometry on ∆ϕ, the TAWSS, and the TWSSG is enhanced when the Womersley number α and/or Reynolds number Re increase. Compared to the straight vessel, the phase shift increases at the inner bend of the curved vessel and inside the aneurysm dome. The strongest effect of geometry on ∆ϕ was observed in a curved tube containing an aneurysm, where the spatial behavior of ∆ϕ becomes highly irregular (Figure 6c). This behavior must be caused by the complex flow field, including nontrivial backflow effects, that prevails in the aneurysmal dome.
Importantly, the geometry considerably affects the connection between ∆ϕ and other WSS-based metrics, such as the TAWSS and the TWSSG (Figures 6 and 7). Given that the imposed signal is a sine wave, it was expected that the TAWSS and TWSSG would provide qualitatively the same type of information (see the similarity of the distribution of the TAWSS and TWSSG in Figure 6).
Besides the strong similarity between the TAWSS and TWSSG in all investigated geometries, we found that ∆ϕ and the TAWSS in the curved tube are closely correlated as well: ∆ϕ is large in places where the TAWSS is high and small or even negative where the TAWSS is low (Figure 7a,b).  (Figure 6c). In (c), large WSS appears in the flow impingement region. Outside the aneurysm, there is a noticeable positive correlation between ∆ϕ and TAWSS. Within the aneurysmal dome, however, no strong correlation is found.
The correlation between ∆ϕ and the TAWSS seems to be lost as the flow inside the aneurysm dome becomes less regular. A possible rationale for this observation could be that velocity fluctuations are enhanced due to the nonlinearity of the fluid dynamics equations in the investigated inertial flow regime. This interpretation is in line with the observation that the phase shift grows as the Reynolds number increases (see Figure 5). Velocity fluctuations may be of a purely numerical nature, but may also arise from the irregularities of the geometric structure. Since the correlation between the shear stress and phase shift persists in the case of the curved tube without aneurysm, we concluded that fluctuations arising from geometry (possibly augmented by grid resolution effects) are the dominant reason for the loss of correlations in the aneurysm dome (Figure 7c).
It seems that ∆ϕ and the TAWSS can be used as proxy metrics for each other, as long as the geometry and flow field are sufficiently smooth. Contrarily, the flow field in an aneurysmal dome is complex, and the TAWSS and ∆ϕ show different behaviors. In principle, the combined effects of a low TAWSS and a large variation of ∆ϕ in the aneurysm may trigger biological processes that could impair the mechanical stability of the endothelium over a long period of time. A more careful examination of the sensitivity of these biological processes seems justified.
Flow is sensitive to geometry, which can, in turn, induce significant differences in the phase shift, ∆ϕ. If ∆ϕ were determined solely by the primary (axial) flow, ∆ϕ would be essentially the same everywhere. Due to the observed dependence of ∆ϕ on Re and tube curvature, it seems that the spatial variation of ∆ϕ is closely related to secondary flows [54] that appear in complex geometries when inertia is not negligible. Due to inertia, a perturbation in pressure takes some finite time to be visible in the surrounding velocity field, and consequently in WSS. Thus, for simple geometries such as a straight tube, WSS lags behind the pressure and ∆ϕ > 0 (Figures 3 and 6a). Yet, ∆ϕ < 0 can occur due to geometrical effects. While laminar flow in a straight tube is perfectly axial, in a curved tube (Figure 6b), the curvature leads to secondary flows within the plane perpendicular to the axial direction. This is illustrated in Figure 8a, which shows the cross-section of the curved tube midway between the inlet and outlet. We observed that the magnitude and features of the secondary flows depend on the tube curvature and the flow properties, in particular Re, which affect the mechanism of momentum transport (data not shown). For flow in a straight tube, WSS is determined purely by the radial diffusion of momentum due to viscosity. When the geometry is curved, however, WSS can also be influenced by the diffusion and advection of momentum via secondary flows.
The curvature of the tube induces centripetal forces in the fluid, which move the line of maximum velocity away from the centerline of the tube towards either the bottom or the top of the tube. From Figure 8a, it can be seen that the flow is shifted towards the bottom of the cross-section, leaving the top wall with a lower velocity and WSS (see the TAWSS in Figure 6b). Figure 8b tells us that, in the presence of curvature, the delay between pressure and velocity signals is the largest for the region within the channel with the maximum velocity (the rightmost data in Figure 8b). Moreover, the variation of ∆ϕ when approaching the vessel wall is nonsymmetric, as it seems to saturate on a plateau towards the bottom wall, but decreases approximately linearly when moving towards the upper vessel boundary. There exist three possible sources of velocity (or WSS) variation: radial diffusion of momentum, radial advection of momentum, and centrifugal acceleration (or force). The last term drives the secondary flow and acts as an additional pressure gradient. We speculated that the phase shift between the pressure and WSS signals is controlled by the combined effects of these three sources. Although the detailed mechanisms underlying the phase shift remain unclear, we can conclude that the axial advective transport of momentum no longer dominates the phase shift, provided that the pipe (or vessel) is not straight.
We found that the observed phase shift is a complicated function of velocity field, geometry, and flow parameters (Re and α). There do not seem to exist simple scaling laws for ∆ϕ that could be used to predict the local phase shift in a given geometry without running a full simulation of the flow field. At the current stage, a predictive model for ∆ϕ is not yet available, and the present work thus emphasizes the need for further research.  [53]. Farther away from the site of highest velocity, the phase shift is smaller in the straight tube; this also applies to the curved pipe in addition to a small region above the site of highest velocity. Note that the point of highest velocity is defined as the point on the cross-section that has the largest velocity averaged over a full period. The peak velocity at any given point indicates the maximum velocity within a full period at that point.

Conclusions
The phase shift (or time delay) between pressure and wall shear stress (WSS) signals can elicit pathological responses in endothelial cells and play a crucial role in diseases, such as atherosclerosis and brain aneurysms. In this work, we employed a hybrid immersedboundary-lattice-Boltzmann method to investigate the phase shift, with the focus on the effects of flow oscillation and inertia (via Womersley and Reynolds numbers) and wall geometry (straight tube, curved tube, and curved tube with an aneurysm). The accuracy of the method was shown through pulsatile Womersley flow tests using a straight rigid pipe.
For pulsatile flow in a straight tube, the phase shift between the pressure and WSS increases with the Womersley number, but is independent of the Reynolds number (or flow speed). Within a curved tube, the phase shift shows a strong dependence on the Reynolds number as well, an effect that is caused by inertial momentum transport from the tube centerline to the walls. Importantly, we found highly heterogeneous distributions of the phase shift within the aneurysmal dome.
While there is no simple way to precisely predict how the phase shift behaves in an even slightly complex geometry, we found that the phase shift is well correlated with the time-averaged wall shear stress (TAWSS) and the temporal wall shear stress gradient (TWSSG) in the curved tube, outside of the aneurysm. Interestingly, the strong correlation between the phase shift and the TAWSS/TWSSG is lost within the aneurysmal dome. Our findings, therefore, suggest that the phase shift could be used as an independent flow metric in segments of complex geometry where the influence of the flow on the blood vessel wall is not reliably characterized by WSS alone. We speculate that geometries with an uncorrelated phase shift and WSS patterns might be particularly relevant to long-term detrimental biological effects that weaken the blood vessel integrity eventually.
A quantitative predictive model for the phase shift as a function of geometry and flow parameters (Reynolds and Womersley numbers) is not yet available. The phase shift seems to depend strongly on the details of the secondary flows that arise in curved geometries at finite inertia. At the current stage, a full flow simulation is necessary to predict the phase shift distribution in a given geometry. More research into the underlying physical mechanisms of the phase shift seems warranted.
In this work, blood vessel walls were assumed to be rigid and the effect of red blood cells were neglected. We expect that the flexibility of the vessel wall and the particulate properties of blood may also play a role in the phase shift between pressure and shear stress signals. A thorough study of this issue provides a challenging task for future work. Acknowledgments: H.W. acknowledges his IMPRS-SurMat scholarship and the financial support by the STKS department at ICAMS, the Ruhr-Universität Bochum. ICAMS acknowledges funding from its industrial sponsors, the state of North-Rhine Westphalia, and the European Commission in the framework of the European Regional Development Fund (ERDF). We acknowledge support by the DFG Open Access Publication Funds of the Ruhr-Universität Bochum.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of the data; in the writing of the manuscript; nor in the decision to publish the results.

Appendix A. Time Scales
The pressure signal propagates at the speed of sound. The shear stress has two transport mechanisms: advection and momentum diffusion. For flow in a straight tube, advection is only along the flow direction and therefore does not contribute to the radial transport near the wall where the flow is tangential and slow compared to the tube center. Therefore, in a straight tube, WSS is determined by momentum diffusion. For flow in a curved vessel, however, there can be significant advection from the centerline to an area close to the wall due to inertia and geometry-induced secondary flow. In that case, the overall transport of shear stress from the centerline to the vessel surface will not just be due to diffusion.
Upon a first view, the phase shift might be associated with five physical time scales: momentum diffusion time t ν ∼ R 2 /ν, inertial time t ρu ∼ ρU/P 0 , oscillation time t ω ∼ 2π/ω = T, sound wave (or acoustic) time t s = L/c s , and advection time t u ∼ L/U ∼ t s /Ma. Here, R and L, respectively, indicate the radius and length of the tube or blood vessel; ν is the kinematic viscosity, ν = η/ρ, where η and ρ, respectively, represent the dynamic viscosity and the fluid density; U indicates the mean flow velocity; P 0 is the amplitude of pulsatile pressure gradient; ω is the angular frequency (ω = 2π/T with T indicating the period); c s is the speed of sound; Ma denotes the Mach number.
For pulsatile flow in a straight rigid pipe, the total mean (i.e., spatiotemporal average) velocity U is U = P R 2 /(8η) where P indicates the time average of the pressure gradient [52], e.g., P = 1 T T 0 P (t) dt = 1 T T 0 P 0 (A sin(ωt) + B) dt = BP 0 where A and B are constants. Thus, the inertial time can be written as t ρu ∼ ρU/P 0 = P R 2 /(8νP 0 ) = BR 2 /(8ν) = t ν B/8, i.e., t ρu ∼ t ν . Hence, in the present problem, there are only four independent times: the momentum diffusion time t ν , the oscillation time t ω , the advection time t u , and the acoustic time t s .
The speed of sound, c s , in blood is roughly 1570 m/s. The radius, R, of major cerebral arteries is about 2 mm. In experiments or simulations, a length-to-diameter ratio no less than five is generally a good choice. In this work, we used L/R = 12. Blood density ρ and dynamic viscosity η are normally within the range of 1-1.06 g/cm 3 and 3-4 mPa s, respectively. The period of the human heartbeat (or the period of pulsatile pressure), T, is 0.8 s in normal cases. Given the above parameters, the traveling time of the pressure signal along the short pipe considered in this work was only 2 · 10 −5 T, which is negligible compared to all other physically relevant times.