Physiologic blood flow is turbulent

Contemporary paradigm of peripheral and intracranial vascular hemodynamics considers physiologic blood flow to be laminar. Transition to turbulence is considered as a driving factor for numerous diseases such as atherosclerosis, stenosis and aneurysm. Recently, turbulent flow patterns were detected in intracranial aneurysm at Reynolds number below 400 both in vitro and in silico. Blood flow is multiharmonic with considerable frequency spectra and its transition to turbulence cannot be characterized by the current transition theory of monoharmonic pulsatile flow. Thus, we decided to explore the origins of such long-standing assumption of physiologic blood flow laminarity. Here, we hypothesize that the inherited dynamics of blood flow in main arteries dictate the existence of turbulence in physiologic conditions. To illustrate our hypothesis, we have used methods and tools from chaos theory, hydrodynamic stability theory and fluid dynamics to explore the existence of turbulence in physiologic blood flow. Our investigation shows that blood flow, both as described by the Navier–Stokes equation and in vivo, exhibits three major characteristics of turbulence. Womersley’s exact solution of the Navier–Stokes equation has been used with the flow waveforms from HaeMod database, to offer reproducible evidence for our findings, as well as evidence from Doppler ultrasound measurements from healthy volunteers who are some of the authors. We evidently show that physiologic blood flow is: (1) sensitive to initial conditions, (2) in global hydrodynamic instability and (3) undergoes kinetic energy cascade of non-Kolmogorov type. We propose a novel modification of the theory of vascular hemodynamics that calls for rethinking the hemodynamic–biologic links that govern physiologic and pathologic processes.

The pioneering work of Womersley 1 and his coworkers 2,3 essentially laid the foundation of modern hemodynamics research. By developing a time-dependent one-dimensional exact solution of the Navier-Stokes equation, Womersley showed that blood flow in main arteries can be described by a Fourier decomposition of the cardiac harmonics 3,4 . This work has been later extended to account for wall elasticity 5,6 and non-Newtonian blood viscosity 7 . The Womersley flow model (WFM) became a founding principle upon which modern blood hemodynamics studies are based. Researchers have assumed, based on WFM, that blood flow is essentially laminar, and transition to turbulence (or presence of disturbed flow, which is a poorly defined hemodynamic term often used in medical research) shifts the blood hemodynamics leading to the initiation of vascular diseases such as brain aneurysms or atherosclerosis 8,9 . Such disease association is attributed to the mechano-sensing properties of endothelial cells that make them responsive to various flow properties 10 . Therefore, the accurate identification of blood flow regimes is an essential step in characterizing the hemodynamic patterns that govern endothelial cells mechanobiology 9 . We have previously shown how the inaccurate assumptions of blood viscosity and the misinterpretation of wall shear stress (WSS) 9,11 have impacted intracranial aneurysm research, leading to inconsistently varying and contradicting results 9 . Moreover, we have recently shown that turbulence exists in pulsatile multiharmonic flow of mean Reynolds number Re m ∼ 300 in idealized model of intracranial aneurysm flow using particle imaging velocimetry (PIV) 12 . Similar transitional and turbulent regimes were detected in intracranial aneurysm by other research groups in vitro 13 and in silico [14][15][16][17][18] . In a seminal article, Jain et al 18 showed that at peak systolic conditions, intracranial aneurysm exhibit random velocity fluctuations and kinetic energy cascade at Re < 400 in the parent artery. Subsequently, Jain et al 17 argued that such complex flow patterns may alter our understanding of aneurysm growth and rupture. Turbulence casts complexity not only on the hemodynamics of intracranial aneurysm, but also on such of carotid occlusive disease 19,20 . It is established that disturbed flow is

Results
Here, we present the comparative evidence on the inherent turbulence in the carotid artery and we argue that similar behavior should be observed in other major arteries of the human body. Results from the exact solution of three other arteries (aortic root, thoracic aorta and iliac) are also presented. We hypothesize that the intrinsic multiharmonic waveforms of blood flow have properties that drive blood flow towards chaotic turbulent state. As much as it is important to characterize the differences between vessels in such regard, the authors believe that it would be cumbersome to conduct such characterization in this article. We invite the community to reproduce the following results using other sets of data for different arteries to further investigate inherent blood flow turbulence.
Vascular blood flow is chaotic. First, we borrowed tools from the chaos theory to test for the SDIC 27 , and to demonstrate that blood flow, expressed by both exact solution and in vivo measurements, is not periodic. We have obtained Lyapunov exponents, as a test of SDIC associated with turbulence 28 , from the time series of blood   Figure 2a shows the trace of Lyapunov exponents in orbital space for the carotid artery. Similar Lyapunov exponents obtained from in vivo DU measurements are shown in Fig. 2b,c for two healthy volunteers (two of the authors). Positive Lyapunov exponents, indicating SDIC, has been previously reported in few studies investigating pulsatile flow in an artificial heart 29 . However, this is the first evidence on the existence of SDIC in Womersley equation. A more elaborative discussion of relevant researches on SDIC of blood flow is presented in the discussion section.
Vascular blood flow is globally unstable. The theory of hydrodynamic stability has been developed to evaluate the stability of steady flow under finite space-time perturbation 30,31 . Some works have evaluated the stability of monoharmonic (i.e. sinusoidal) flow [32][33][34][35][36] . However, there are no universal criteria that can generally describe the stability of multiharmonic pulsatile flow. We have extended the available criteria to evaluate Womersley flow. In Womersley flow, the velocity field is given as: where n is the total number of harmonics in a specific waveform. Here, we propose that a given Womersley flow is globally unstable if lim Figure 3 shows colourmap of such condition as obtained both from exact solution and in vivo DU measurements of the carotid artery. It is clear that the condition is achieved in both datasets confirming that blood flow is globally unstable. In Fig. 4, the acceleration field obtained from the exact solution is depicted. Maximum acceleration and minimum deceleration are proximal in time domain, while the remaining time domain of the acceleration field is dominated by quasi-steady flow. This supports the assumption made to evaluate the global instability, bearing in mind that Physiologic vascular blood flow exhibits kinetic energy cascade. One of the most profound properties of turbulence is the cascade of kinetic energy. Kinetic energy is transferred from larger to smaller vortex structures and ends up irreversibly dissipating to the surrounding media in the form of heat. Figure 5 shows kinetic energy cascade in carotid blood flow as depicted from exact solution and in vivo DU measurements. To the best of the authors' knowledge, the depicted slopes of such cascade do not match any of the Kolmogorov scales reported in literature. Turbulence, in such case, is of non-Kolmogorov regime. Such uncharacteristic turbulence regime has been observed in atmospheric turbulence 37,38 , however, it is the first time to be reported in this physiologic setting. Until the present moment, the main mechanism linking hemodynamics to biological and pathological processes in arteries is the wall shear stress 8 . The existence of kinetic energy cascade in blood www.nature.com/scientificreports/ flow should develop this mechanism to include other mechanical and kinetic factors that could contribute to mechanosensory and mechanotransduction of endothelial cells.

Discussion
In fluid mechanics, laminar flow is a special case in which the nonlinearity of Navier-Stokes equation is far less than such of the general case of turbulent flow. It is difficult to trace the origins of blood flow laminarity assumption in the early works of Womersley et al. However, it is generally believed that since the WFM describes quasi-periodic flow in cylindrical conduits (i.e. vessels) with Reynolds number below 1,000, it would not fit into the Kolmogorov-Obukhov statistical theory of turbulence. The Reynolds criteria was used to classify vascular blood flow, although its essential assumptions (i.e. steadiness, Newtonian viscosity and uniform geometry) do not apply on blood flow. Hence, the classification of Womersley flow was established to be laminar flow in normal human arteries 39,40 . In the following decades, all the vascular diseases models were established based on that classification. The deviation from the so called laminar flow conditions is believed to be linked to atherogenesis and aneurysm formation through various mechanisms associated with endothelial cells inflammatory reactions and dysfunction 8,[41][42][43][44] . This reasoning has dominated vascular biology research, and most of the body of literature regarding endothelial cells response to hemodynamics for the past four decades 8,45 .
The notion of turbulence in relevant literature often refers to Reynold's turbulent flow criteria 46,47 that describes fully developed turbulent flow. The associated interpretation of turbulence is, therefore, limited to isotropic homogenous turbulence. Kolmogorov 48,49 and Obukhov 50,51 proposed the statistical theory of isotropic homogenous turbulence based on fully developed turbulent flow 52 . Such theory is derived based on key assumptions. Taylor's frozen turbulence hypothesis 53 is one example of these assumptions. It relates spatial statistics to temporal  www.nature.com/scientificreports/ statistics and holds only if ũ/U ≪ 1 . In our work, such ratio has larger ranges, as shown in Supplementary Figure S1. Therefore, the notion of turbulence in our work indicates turbulent flow that does not fall within the Kolmogorov-Obukhov theory. Subsequently, we have borrowed the expression non-Kolmogorov turbulence 54,55 from atmospheric turbulence literature to describe our observations. Despite the differences between the former and the latter, both phenomena cannot be described by Kolmogorov-Obukhov theory. Turbulent flow has many characteristics that keeps it as one of the standing mysteries in classical physics 56 . Some of these characteristics are chaos, instability, and kinetic energy cascade 57 . Turbulent flow is sensitive to initial conditions, which makes it subjected to the properties of chaos theory 58 . A turbulent field is instable where any finite perturbation in space and time propagates. The vorticity field is always non-zero in turbulent flows, and the kinetic energy transfers from large to small vortices and vice versa 59 . The cascade of kinetic energy is a characteristic feature of turbulent flows and it contributes to the very definition of turbulence 60,61 . Steady, periodic and quasi-periodic laminar flows tend to conserve their kinetic energy due to the absence of the vortex formation and breakup phenomena. The classical Kolmogorov-Obukhov theory of turbulence describes isotropic homogenous turbulence in which energy cascade is subjected to scaling laws that define the rate of transfer of kinetic energy 62 . Figure 6 shows the kinetic energy cascade, in frequency domain, for four arteries at two different radial locations based on the WFM exact solution. It is clear from Fig. 6 that the slope of kinetic energy cascade is affected by spatial location (i.e. radial distance from the wall). This is further explained in Supplementary Figure S2, where the radial profiles of kinetic energy cascade slope are plotted. Supplementary Figure S2 shows that the slopes vary with radius, in the four arteries, with values close to Kolmogorov scaling laws (  www.nature.com/scientificreports/ in vivo data. Visee et al. 71 used consecutive transcranial Doppler (TCD) to investigate the existence of chaos in patients with occlusive cerebrovascular disease. They had evidently shown that blood flow in healthy conditions exhibited positive Lyapunov exponents, while in impaired intracranial circulation it was found to be of periodic nature with near-zero Lyapunov exponents. These findings have been also supported by the works of Ozturk 72 and Ozturk and Arslan 73 where TCD signals showed positive Lyapunov exponents in the intracranial circulation. In our analysis of the Womersley equation and in vivo DU measurements, we performed Lyapunov exponent analysis, which essentially informs us whether the data or equation are inherently chaotic or not. Womersley equation had a positive Lyapunov exponent, meaning that it has properties of chaos under multiharmonic pulsatile flow that is the blood flow. This, by definition, negates the assumption that the WFM describes a pulsatile laminar blood flow regime. Thus, we argue that a transition to turbulence per se is not a cause for disease initiation, rather a change in turbulence characteristics due to complex geometry or harmonics would be the culprit. This of course adds a layer of complexity to vascular hemodynamics and biology. In fluid mechanics, it is generally believed that hydrodynamic instability is always associated with turbulence 74 . However, most of the work done to understand the correlation between the two phenomena in pipe flow was mostly limited to steady flows subjected to finite perturbations in space and time 75 . The stability of pulsatile flows have been studied in some works [76][77][78][79] , however, there is no consensual methodology among fluid dynamicists for investigating their stability. Therefore, to study the global instability of blood flow both from Womersley equation and in vivo measurements, we had to adopt the main principles of the hydrodynamic stability theory and develop intuitive representation of its criteria for global instability. The concept underlying this approach is fairly simple. Serrin 80 argued that for a viscous flow to be considered stable "the energy of any disturbance tends to zero as t increases". Hence, we assumed that each harmonic component of the flow waveform (n) could be viewed as a space-time perturbation to its predecessor in in frequency domain (n − 1) , as shown in Fig. 3. In Fig. 7, the stability condition derived from this principle is plotted for four arteries based on the exact solution of Womersley equation using boundary conditions from HaeMod database. We have used 40 harmonics to represent the solution of Womersley equation. The reasoning behind this is based on a mass transfer analysis. We have analyzed the contribution of each harmonic in mass transfer by computing the mass flow rate, per unit density, ṁ = A × u(x i , t) with a loop that changes the number of harmonics (n) from 1 to 1,000. We have found that lim n→∞ṁ n+1 −ṁ ṅ m n = 0 . Supplementary Figure S3 shows this difference for the first 40 harmonics. Based on When one thinks of laminar pulsatile flow, the only available flow field variable relevant to physiological vessel changes is wall shear stress, which is well studied and characterized especially in in vitro endothelial flow exposure experimentation 10 . On the other hand, when turbulent flow is considered, numerous variables could affect the vessel on multiple scales. For instance, kinetic energy, represented by the energy cascade, transfers from large coherent vortices to small isotropic homogenous vortices before dissipating in viscosity scales as heat. We have evaluated the length scales and energy scales of such vortices and they fall well within the limits of cellular mechano-sensing 81 . This begets an important question, how do different energy cascades, for example a direct versus an inverse cascade, impact endothelial cells? and more importantly, how can we fine tune our experimental setups to measure and reproduce these conditions? Moreover, how would these new findings impact vascular hemodynamics based therapeutic such as shear driven drug delivery, a possible drug delivery method to areas of altered shear stress 82 , or intravascular devices 83 ? All of these questions, indeed, are interesting prospects of the presented results.
The main obstacle in trying to answer these questions from the available literature, is that up to this point, turbulent flow (referred to in most of the literature as disturbed flow) is considered to be pathologic as opposed to laminar flow 84,85 . Moreover, on most studies that discussed the differences between laminar and turbulent flow regimes on endothelial mechanobiology, both laminar and turbulent flows were at different WSS values 8,86,87 . Given the fact that WSS is one characteristic of flow, and the complexity of turbulent regimes, making a change of WSS impactful within various aspects of the turbulent flow that might in turn impact endothelial mechanosensing via unidentified mechanisms, it is difficult to draw direct comparisons between the impact of laminar and turbulent flows starting from the assumption that turbulent flow itself is physiologic in healthy carotid arteries in absence of stenosis as we propose herein. To further elaborate on this notion, we previously performed a study on endothelial cells using laminar flow only at various WSS values and analyzed the microRNA expression using microarray 10 . To our surprise, we observed new unreported microRNAs to be differentially regulated between high and low WSS conditions, while previously reported microRNAs were not significant in our analysis. When we consider that most literature studying endothelial microRNA in flow responses uses turbulent flow at lower WSS compared with laminar flow at physiologic WSS 45 , the importance of performing parametric studies to characterize the differences between turbulent and laminar flow at same or similar WSS becomes obvious. www.nature.com/scientificreports/ In conclusion, this work aims to propose a modification of the theory of vascular hemodynamics. That is, the blood flow is inherently turbulent and not laminar, and changes in turbulence kinetic energy or other properties of turbulence could be the driving factors behind the hemodynamics-biology links. These results and their theoretical consequences should motivate the scientific community to update their thinking regarding hemodynamic drivers of endothelial and vascular processes, given the inherent complexities and chaos associated with turbulence. Moreover, more hemodynamic research, with accurate and rigorous methodologies, should aim at further characterizing the interesting features of turbulent blood flow in physiology and pathology. To improve the accuracy of WFM, similar analysis should also be conducted with mathematical and computational models that accounts for non-Newtonian effects, platelet and blood cells dynamics and fluid-solid interaction with arterial wall. where k s e iωt is the oscillating pressure gradient, r is the artery radius, µ is the viscosity of blood, � = r ρω µ is the Womersley number, J 0 is the Bessel function of zero order and first kind,

Materials and methods
is the complex frequency parameter, and ζ (r) = � r R is the complex variable. The derivation of Eq. (2) from Navier-Stokes equation can be found in the original paper by Womersley 4 or and in more instructive details in 88 . The exact solution of Eq. (2) as initial-boundary value problem has been extensively reported in literature, however, with no evidence of checking the existence of turbulence and chaos. In this work, the initial-boundary conditions were obtained from the open-access database HaeMod (https :// haemo d.uk/virtu al-datab ase). This database provides blood flow waveforms for 11 major arteries. The present work establishes the investigation on the carotid artery (using the database as well as via doppler U/S measurements from volunteers) and provides the exact solution results for three more arteries in supplementary materials. The waveforms were created using physiologically realistic model based on human data. HaeMod database has been extensively validated in literature and confirmed to represent physiologically relevant and valuable data. However, Fourier decomposition revealed that the waveforms are limited to 40 harmonics. No institutional approval was required for the DU measurements as per Tohoku University and Kohnan hospital (the location of the DU measurements) ethics guidelines.

Lyapunov exponent calculations.
For any time series dataset, the rate of separation of infinitesimally close orbital trajectories are characterized by Lyapunov exponent where the initial separation is δX o and the divergence rate is |δX(t)| ≈ e t |δX o | where t is time and is the Lyapunov exponent. We have used the opensource code provided by Wolf et al. 89 and described in their 1985 famous paper 27 to build the phase-space and orbits of the time series datasets obtained from the exact solution of the Womersley equation and DU measurements.
Evaluation of global hydrodynamic stability in pulsatile flow. The onset of transition in physiological flow requires hydrodynamic instability to commence. Such instability is caused by a sustained disturbance in the flow field, where its energy is expressed as: is the blood velocity waveform with a i and b i as the Fourier coefficients of the Womersley waveform, and u o (x i ) is the steady component of the flow which is characterized by Hagen-Poiseuille parabolic velocity profile. The sustenance of a disturbance must be global (i.e. in space and time) for transition to take place.
In order for the flow to become globally instable, the time rate of change of E V must be positive, hence the disturbance can possibly prevail over the viscous resistance of the fluid. The Reynolds-Orr equation 90 for disturbance energy reads: www.nature.com/scientificreports/ where t is time, S ij and s ij are the strain rates of basic and altered (i.e. disturbed) flows, respectively, Re is the Reynolds number and V is the wall-bounded flow volume. The first term on the right-hand side (RHS) of (3) describes the exchange of energy with the basic flow, while the second term describes the energy dissipation due to viscosity. If the former is higher than the latter (i.e. the disturbance energy does not decay with time), the flow will be globally instable. Thus, the criteria for global stability can be written in terms of critical Reynolds number as: In Eq. (5), the critical Reynolds number Re cr,m defines the global and monotonic stability criterion for viscous fluid flow subjected to a disturbance of the velocity u(x, t) . Hence, if Re > Re cr,m the flow becomes monotonically instable, in other words the disturbance will not decay for t → ∞ and transition will take place. Hence, asymptotic instability occurs when lim n→∞ E V (t) E V (0) � = 0 . In pulsatile multiharmonic flow, such condition takes the form lim n→∞ E V (n) E V (n−1) � = 0 where n is the index of harmonic in the Womersley velocity waveform described in Eq. (3). This condition has been evaluated in space and frequency domain using Eq. (3) and the trapezoidal method for numerical differentiation. Following Thomas 91,92 and Hopf 93 , Serrin 80 proved that in order to have monotonic instability in any viscous flow, the disturbance energy must satisfy the following condition: where E VO is the initial disturbance kinetic energy, U mc is the velocity of the basic flow during a one cycle of the pulsatile flow. Serrin's second theorem proved that Re cr,m = 32.6 for viscous incompressible wall bounded flows. Reynolds number in main arteries is larger than the latter value. Hence, blood flow, by definition, is monotonically instable.
Calculating energy cascade in frequency domain. The local kinetic energy in frequency domain was calculated as 94 : where f is the frequency, F is the fast Fourier transform and L the length of ũ i matrix. It was not possible to use Taylor's hypothesis to calculate the kinetic energy in wavenumber domain www.nature.com/scientificreports/ since ũ U > 1 . Figure 8 shows velocity patterns in radial and time domains of the carotid artery to show that the conditions for using Taylor's hypothesis do not apply.

Data availability
Raw data required to reproduce the carotid artery analysis based on WFM has been deposited at Harvard Dataverse and is available for public access under CC license at: https ://doi.org/10.7910/DVN/OM1T4 2 Received: 17 April 2020; Accepted: 27 August 2020