Theory to Predict Shear Stress on Cells in Turbulent Blood Flow

Shear stress on blood cells and platelets transported in a turbulent flow dictates the fate and biological activity of these cells. We present a theoretical link between energy dissipation in turbulent flows to the shear stress that cells experience and show that for the case of physiological turbulent blood flow: (a) the Newtonian assumption is valid, (b) turbulent eddies are universal for the most complex of blood flow problems, and (c) shear stress distribution on turbulent blood flows is possibly universal. Further we resolve a long standing inconsistency in hemolysis between laminar and turbulent flow using the theoretical framework. This work demonstrates that energy dissipation as opposed to bulk shear stress in laminar or turbulent blood flow dictates local mechanical environment of blood cells and platelets universally.


Introduction
Turbulence is ubiquitous, and is particularly prominent in engineered cardiovascular devices as well as pathophysiological blood flow. There is strong evidence that turbulence impacts the environments of erythrocytes and platelets on the cellular level [1,2] in a manner physically distinct from that known in simple laminar shear flows such as in a viscometer. While our physical understanding of the structure of turbulence and its universal properties has received significant attention in the past [3][4][5][6][7], there is no theory that links the statistical properties of turbulence to shear stresses physically experienced by cells transported in whole blood flow. Shear stress acting on cell membranes is a critical mechanical cue that regulates biological activity [7][8][9][10] and therefore a theory relating turbulence to shear stress environment of cells is necessary.
In turbulent blood flow, the complex spatio-temporal fluctuations of shear stress leads to hemolysis and platelet activation [11,12]. Such phenomena are critical when designing life saving devices such as artificial hearts, ventricular assist devices, stents, grafts, and heart valves. The phenomena can further impact disease such as in the case of an aortic stenosis or atherosclerosis. Current models that predict stress experienced by blood cells are purely empirical and based on classic experiments [12] that yield paradoxical results under laminar and turbulent conditions [13]. For a comprehensive review of studies on turbulent blood flow related hemolysis and platelet activation, the reader is directed to Refs [12,[14][15][16][17]. As discussed in Ref [12], despite many of these studies, there has not been a solid physically justifiable connection made between turbulence and the shear stress acting on blood cells and platelets. A strong requirement of a physical theory is that it should make a link, in a manner independent of laminar and turbulent regimes of flow because the pertinent parameter is flow at the length scale of individual cells, where it is considered laminar. Another aspect that is important to consider is the notion of universality of turbulent structures despite the intermittency issue [6]. Do the universality properties of turbulence hold in the most complex of blood flows considering Newtonian and non-Newtonian properties of blood? If so, will the distribution of shear stress acting on blood cells and platelets be universal? Here the term ''universal'' is used not to imply homogeneous isotropic turbulence (HIT), but rather loosely to emphasize the significant agreement between the distributions of instantaneous dissipative scales in complex in-homogeneous shear flows with the distributions observed in HIT as well as the robust presence of inertial range scaling [18][19][20][21].
The goal of this work is to address the above by introducing a unified (in the sense of laminar vs. turbulent) physical theory to: (1) identify the relevant dynamic properties of flow that link to the predicted shear stress experience by cells based on fundamental physical arguments. We achieve this for the special case of blood, but the underlying physical arguments may hold for any cell suspension constrained to the Newtonian regime of its rheology, (2) theoretically consider the relevance of non-Newtonian effects on the smallest scale turbulent structures for blood based on order of magnitude arguments; and (3) test the ''universality'' of small scale structure of turbulent flows in one of the most complex turbulent blood flow problems (flow through a bileaflet mechanical heart valve) and experimentally examine the universality of the predicted shear stress distribution, thus testing the robustness of the theory. Since we lack the capability to experimentally measure shear stress on individual cells transported in a turbulent flow, we use the new theoretical framework to resolve inconsistency in published hemolysis data in laminar and turbulent pipe flow as a surrogate validation of the theory. Finally, we note that the theoretical framework presented in this work is not limited only to blood cells, but applies to any case of suspended cells that are sufficiently smaller than the smallest scales of turbulent motion.

Methods: Theoretical Construction
Theoretical development is initially constrained to turbulent flows where the smallest turbulent eddies are greater than the size of the red blood cell (RBC), i.e. the instantaneous minimum size of turbulent eddy is always .O(10 mm). This is an important limit, as it is not possible to physically represent turbulent eddies smaller than the size of the biggest cells with the single-phase approximation of blood. In fact, the single-phase approximation may very well be susceptible to errors even when local eddies are within two or three times the size of the red-blood cell, i.e. ,25 mm. For now, let us accept this as a limitation and later discuss repercussions of this assumption. Nevertheless, for turbulent eddies .25 mm, the single-phase continuum representation of blood may be considered valid [12]. It is also important to underscore that using a continuum single-phase model of blood cannot be equated to the assumption of RBCs passively seeded on to a continuous medium.
On the contrary, the single-phase representation of blood ''lumps'' all effects of the physical reality into the constitutive rheology of blood.
The theoretical construction is next focused on whether small scale (dissipative) turbulent structures may be considered Newtonian or non-Newtonian. This is done through order of magnitude arguments combined with existing knowledge about blood rheology. Subsequently, we introduce basic arguments with respect to the nature of turbulent scales of motion in turbulent blood flows expected in the cardiovascular system and artificial devices. Finally, physical arguments are introduced that link turbulence statistics to the distribution of laminar shear stresses acting directly on cells.

Newtonian vs Non-Newtonian
It is well known that whole blood significantly deviates from Newtonian behavior when local bulk shear rates are below the order of 100 s 21 [22][23][24] or when a vessel diameter is less than 100 mm [25]. Now, consider the instantaneous turbulent dissipative length scale, g, defined such that the local instantaneous velocity increment (difference), du g , across two points in space separated by g is such that the local Reynolds number defined as  [18][19][20]24,26,27]. Note that no assumptions of local isotropy or homogeneous turbulence are being made in this framework. For such a small dissipative length scale, which in a way is the true measure of the actual size of the smallest eddy, the shear rate within such an eddy can be estimated as _ c c g * du g g * n g 2 . Thus the shear rate within such an eddy is set purely by n and g. Now, non-Newtonian behavior of blood dictates that n is dependent on _ c c g , which is physiologically O(1) cSt. It is easy to show that the range of _ c c g * n g 2 corresponding to eddies of Based on these order of magnitude estimates, it follows that instantaneous turbulent eddies in the range 10mmvgv100mm in any cardiovascular blood flow will always be in the Newtonian regime considering _ c c g wO 100s {1 À Á . In other words, any consequence or damage occurring to the cells within these eddies happens while the eddy itself behaves as a Newtonian fluid.

Scales of Motion
With the notion of blood being Newtonian for turbulent eddies 10 mm and above, let us examine the scales of turbulent motion for the smallest scales, without the assumption of isotropy or homogeneity. Recall that turbulent flow is characterized by complex spatio-temporal fluctuations, u i , superimposed with the mean velocity field, U i . These spatio-temporal fluctuations have been phenomenologically represented as turbulent eddies of varying length scales. It is important to clarify that physically there are no circular ''eddies'', but the term eddies only refer to the Fourier analogy of representing the fluctuating field as an summation of kinetic energy over ''eddies'' ranging from the smallest possible scales to the integral length scale. The turbulent kinetic energy (per unit mass), k~1 2 Su i u i T, is transferred across from large to the small scales (so called energy cascade) with viscous dissipation occurring at the smallest eddies. The Kolmogorov length scale, g K : n 3 =SET À Á 1=4 , where n is the kinematic viscosity and e is the average local kinetic energy dissipation rate (per unit mass), corresponds to the characteristic eddy size based on e. The velocity scale of this eddy equals u g K~n =g K , and therefore the Reynolds number, Kolmogorov eddy is unity. In the context of blood damage, several studies have focused on the relevance of the Kolmogorov eddy to predict the capacity of turbulent flow to damage blood cells (reviewed in [12]). While g K in theory represents the average size of the dissipative eddies (for the case of HIT), it does not, in itself, reflect the smallest eddy in a turbulent flow in realistic turbulent flows. In fact, there exist eddies that are a fraction of the Kolmogorov eddy, so called sub-Kolmogorov eddies, arising from the inherently intermittent nature of the instantaneous energy dissipation rate field [5,6]. Intermittency is a property of all turbulent flows which broadly refers to the fluctuations of the fluctuating velocity gradient tensor s ij : consequently for e~2ns ij s ij . With energy dissipation rates proportional to the square of velocity gradient tensor, the fluctuations in e are far greater and intense, with very large departures from its average. Physically these large positive fluctuations in energy dissipation occur due to local vortex stretching that generate momentarily intense shear regions over very small length scales, much smaller than the Kolmogorov scale itself. In other words, when the energy dissipation rate is locally above the mean energy dissipation magnitude, the corresponding eddy size is sub-Kolmogorov. Thus, the true depiction of micro-environment of cells in realistic turbulent flows corresponds to their tumultuous experiences in a spectrum of dissipative eddies ranging from sub-Kolmogorov scales all the way to the Taylor micro-scale that demarks the largest of the dissipative eddies [6]. Recent turbulence literature estimate the smallest sub-Kolmogorov scales to be of roughly Re 1=4 L times smaller than the Kolmogorov scale [23,26,28]. Here Re L~k 1=2 L n is the local turbulence Reynolds number defined by the integral length scale, L, and the local turbulent kinetic energy. To date, studies addressing turbulence in blood flow have not taken into consideration the issue of intermittency or the sub-Kolmogorov fluctuations thereby ignoring the most damaging aspects of turbulent flow. In this study, we will consider sub-Kolmogorov fluctuations using direct measurements through the definition of the instantaneous turbulent dissipative length scale, g [18][19][20]24,26,27], while constraining our analysis only to flows where the smallest sub-Kolmogorov gw10mm.

Local Energy Dissipation and Shear Stress Acting on
Cells. Based on the above picture of scale distributions without invoking assumptions of isotropy or homogeneity, and introducing the notion of sub-Kolmogorov scales, it is now possible to consider a link between macroscopic properties of blood dynamics or eddy scales to the local shear stress acting on blood cells. Figure 1 illustrates a schematic of RBCs in a hypothetical dissipative scale eddy of size g. Recall that the velocity scale of this eddy is du g~n g to yield the condition of locally unit Reynolds number. Given the length scale and velocity scale, the instantaneous energy dissipation rate corresponding to this smallest eddy is du 3 g g~n 3 g 4 . Looking at Figure 1, let us now assume that almost all of this energy dissipation physically occurs through viscous straining of plasma fluid between cells and that the rate of energy lost in lysing or activating cells is negligible in comparison to heat generation within plasma. This is a valid assumption if (a) it is known that a very small fraction of cells lyse/activate per unit time, and (b) the strain rate of cell membrane deformation is much smaller (by at least an order of magnitude) than the strain rate of plasma surrounding the cell. The relative difference in the strain rates of the cell membrane and the surrounding plasma represents the efficiency of energy transfer from the fluid to strain-energy stored into the cell membrane. Both (a) and (b) seem valid as hemolysis/ activation in most device and physiological flows occur over prolonged duration and the magnitude of free hemoglobin released is relatively very small compared to the total hemoglobin. In particular, the second condition is valid for cases when blood cells are lysed over a prolonged exposure to shearing forces as opposed to instantaneous lysis which has been observed for whole blood shear stresses .450 Pa [29,30]. Complex cardiovascular devices such as the bileaflet mechanical heart valves imposes shear Theory to Predict Shear Stress on Cells in Turbulent Blood Flow PLOS ONE | www.plosone.org stresses in the range ,15 Pa [14]. Furthermore, (b) is supported by the relatively low energy dissipation for both cytoplasm and membrane deformation in a RBC tank treading in flow, corresponding to O(10 215 -10 213 watts) for shear rates O(10 2 ), which extrapolates to O(10 28 watt) for shear rates reaching 10 5 s 21 [31,32].
Given these assumptions, which appear reasonably valid, the rate of energy dissipated within the eddy shown in Figure 1, e, may be approximated to the rate of energy dissipated in the fluid that surrounds cells. The following equation represents the energy balance in watts: Where r B is the density of whole blood, V g is the volume of the eddy, m p is the dynamic viscosity of plasma, _ c c p is the strain rate in the plasma surrounding blood cells, and 0vHv1 is the fractional hematocrit. It is easy to verify that both sides of the equation have units of power. In the above equation, if the plasma dynamic viscosity is known, then it is straight forward to estimate the viscous shear stress, t p~mp _ c c p , within the plasma surrounding the cells as a function of energy dissipation rate. Strictly, the cells experience t p , not t B which can be estimated as t B~mB _ c c B~mB n B g 2 . Clearly, it is evident from a physical basis that the relevant dynamical property of blood flow that dictates t p is the local energy dissipation rate, e.
Recognizing that it is easier to measure grather than e, we solve equation 1 for t p~mp _ c c p after substituting e~n 3 B g 4 to give: The above equation relates the instantaneous dissipative scale g in a general turbulent blood flow to the corresponding shear stress acting on cells. It is important to note that t p and t B are off by a factor t p =t B~m variations in hematocrit under turbulent straining may produce large discrepancies between t p and t B . Figure 2 presents the theoretical estimate of shear stress t p using equation 2 over a range of turbulent eddy sizes 10mmvgv10 3 mm and hematocrit 0vHv0:8. As shown in this figure, t p ranges 10 21 -10 5 dyne/cm 2 . We must note here that we extrapolate equation 1 for g as low as half the size of the RBC. It is interesting to note that for an g about 5-6 mm, equation 1 predicts a shear stress between 400-500 N/m 2 , a magnitude well known to instantly lyse RBCs [29,30]. This illustrates that the t p prediction of equation 1 up g approaching 10 mm appears to asymptotically reach the instantaneous lysis value of 400-500 N/m 2 [29,30]. Furthermore these shear stresses correspond to the release of serotonin from platelets, demonstrating granule release, a process involved in platelet activation [33].
In real turbulent flows, g is a dynamically fluctuating quantity in space and time, around the statistical measure g K : Substituting g K in equation 2 will only represent an average shear  stress acting on RBCs and is perhaps insufficient to reflect the intense fluctuations of the dissipation rate that would correspond to sub-Kolmogorov scale eddies. Thus, to fully capture the statistical nature of shear stress acting on RBC membranes, it is necessary to characterize the probability density function of t p denoted P(t p ). P(t p ) is related to the probability density function of g, Q g ð Þ, through conservation of probability as: may be universal despite the highly intermittent fluctuations of the instantaneous dissipation rate field [18,19]. Analytical forms for this universal distribution exist with good agreement for experiments and simulations [24,27]. The key point relevant here is that P(t p )may be universal in the strongly turbulent regions of blood flow through devices, and estimated using a simple change of variables as: Where The above arguments and physical construction not only provides a mathematical relationship between g(as a surrogate for e), and t p , but also yields P(t p ). However, how does this relationship link turbulence statistics to t p ? The most common and perhaps paradoxical turbulence statistic in the context of blood damage has been the Reynolds stress tensor rSu i u j T. This has been extensively discussed in the past (and key issues reviewed in [12]) with the note that while it does appear to relate to predicting cell damage, there is no consensus with respect to the physical relationship between rSu i u j T and t p . An alternate model to estimate t p based on cell relative velocities between adjacent eddies has been proposed [12] but in our opinion this lacks the physical justification for existence of such intense velocity jumps (more severe than the smallest sub-Kolmogorov event) over submicron scales. Nevertheless, based on the phenomenological arguments made earlier with the hypothetical eddy, it is easy to see that it is the total energy dissipation rate that directly dictates t p . One feasible explanation we can offer to explain the robustness of Reynolds shear stress in predicting blood damage is that for regions in equilibrium, the average rate of energy dissipation, e is related to Reynolds stress given by the equation (proof in [34]): is the mean strain rate tensor.
Plugging this approximation in equation 1 and ensemble averaging it, the root mean square of t p is given by: The above discussion is presented for completeness and offers to reconcile the rather paradoxical issue (so far) of why Reynolds shear stress has been effective in capturing blood damage potential of turbulent blood flows [12]. The reconciliation that we offer through equation 4 is that it is not the Reynolds stress itself, but the product of the Reynolds stress and the local strain rate that determines the energy passed on to the cascade and hence the total energy dissipation rate, which as we have shown should ultimately set the viscous shear stress acting on blood cell membranes. A rough order of magnitude verification of equation 4 can be made based on the contour scale bars in Ref. [14] that experimentally quantifies the Reynolds shear stresses in a heart valve flow. By setting Reynolds shear stress in equation 4,100 N/m 2 (i.e. O(100) N/m 2 ), and substituting m p to be ,10 23 Ns/m 2 and mean strain rate to be O(1000) s 21 , we get St 2 p T 1=2 14N=m 2 . This strikingly agrees with the range of viscous shear stress reported in Ref. [14] to be 15 N/m 2 . This agreement demonstrates further that Reynolds shear stress when combined with the mean strain rate can relate as an order of magnitude estimate to viscous shear stress acting on blood cells. The most comprehensive representation of shear stress acting on cells however is undoubtedly the hypothesized universal distribution function P(t p ).
A Test of Universality and Quantification of P(t p ). With the reasonableness of the Newtonian assumption for physiological turbulent eddies, we interrogate the most complex turbulent blood flow problem -the pulsatile turbulent field surrounding a bileaflet mechanical heart valve using high-resolution phase-locked particle image velocimetry to calculate P(t p ) and assess its universality. Turbulence in mechanical heart valve flows has received enormous attention in medical research with unresolved issues in relation to the applicability of Reynolds shear stresses on blood cells [12]. Briefly, a bileaflet mechanical heart valve was experimentally subjected to physiological aortic conditions and the instantaneous turbulent velocity field was captured in the vicinity of the valve during representative phases of the cardiac cycle. Points of interest within the measurement region where complexity is expected ( Figure 3) were further interrogated to quantify the probability density function Q g ð Þ and the corresponding probability density function of the shear stress,P(t p ) for a hematocrit of 48%.

Particle Image Velocimetry (PIV) of Bileaflet Mechanical
Heart Valve. Experiments are similar to that described in Ref. [35]. The fluid was seeded with polyamide tracer particles (Dantec Dynamics, Denmark) distributed in the 5 to 30 microns range. A laser sheet illuminated the central plane normal to the b-datum line. The laser was generated using the Photonics Industries DM40-527 diode-pump Q-switched laser (Photonics, Bohemia, , NY) with optics to covert the output beam into an expanded laser sheet. The laser had an initial thickness of approximately 1 mm, which was focused down to less than 200 microns within the measurement region using a spherical lens ( f = 1 m). The valve was oriented such that the measurement plane bisected both leaflets at the central plane of valve model. A Photron Fastcam SA3 CCD high speed video camera (Photron, San Diego, CA) synchronized to the laser system via a high speed controller (HSC) (LaVision, Ypsilanti, MI) captured focused images of the illuminated polyamide particles within the laser sheet in the measurement plane. The image area of interest was 1.5D wide and 1D tall with the valve body centered. Image distortion due to curvature of the acrylic tube was compensated in-situ with a calibration plate consisting markers placed in a regular square grid with 1 mm spacing. The DaVis calibration algorithm (LaVision Inc, Germany) automatically tracks the markers and a map to evaluate the corrected image. Corrected image generated of the calibration plate verified successful calibration and distortion correction. The PIV setup achieved a raw data spatial resolution of roughly 27 mm/pixel. PIV measurements were performed in double-frame mode with a laser pulse separation time, Dt = 500 ms. This ensured adequate particle displacements in the range of 10-15 pixels thus maximizing the accuracy of instantaneous velocity measurements to within 2% error. Images were pre-conditioned by first subtracting the minimum image from every image acquired. Instantaneous two-dimensional velocity field was calculated from the raw particle images using cross-correlation processing with a multi-pass scheme. The initial interrogation window size for the multi-pass scheme was at 32632 pixels, which progressively reduced to 868 pixels. Interrogation window overlap was fixed at 50%. Post-processing of the vector data included a median filter that rejected vectors outside 3 standard deviations of the neighbor vector. Gaussian smoothing was used to reduce noise in the vector data. An in-house Matlab code was developed to post-process these raw velocity measurements to derive statistical properties, specifically PDFs.
Peak locking index was calculated to be between 0.02-0.19.
Peak locking index is defined as 4Dw{ 1 4 D where w is the first moment of the probability distribution function of the absolute fractional distance in pixels to the nearest integer pixel displacement. If the probability distribution is uniform in the pixel displacement range 0 to 0.5, then the pixel locking index is zero, indicating no pixel locking. A value of 1 indicates 100% pixel locking (i.e. no sub-pixel displacements recorded). Our range of 0.02-0.19 for peak locking index is far below 0.25, which is the threshold for minor peak locking artifacts. Error Considerations. The sources of error in our measurements are due to resolution as well as random error. Random errors are addressed in this study through statistical averaging of repeated measurements (N = 500) and statistical comparisons. This section briefly outlines the errors in accuracy due to limited resolution of the measurement techniques at hand. The conservative error estimate in velocity is ,2% (i.e. particle displacements may be off by 60.2 pixels out of total displacement of 10 pixels). Laser pulse timing errors are negligible in comparison.
Validation. In order to check if the valve chamber and the acrylic model valve provide equivalent results to clinical quality St. Jude Medical valve, Figure 5 of Ref. [35] compares nondimensionalized leaflet kinematics, and the downstream velocity profile at x/D = 0.33 during peak systole to published results for a clinical quality St. Jude Medical valve [36].
Calculation of Q g ð Þ and P(t p ). PDFs of the local dissipation scale, Q g ð Þ was calculated for each interrogation point in Figure 3. The approach is similar to that previously described [20,37] with the qualifying condition for g being 0:9v gdu g n v2. Briefly, the instantaneous g is estimated from the local velocity increment between PIV velocity measurement points. We have shown in Ref. [37] that a significant portion of the local dissipative scale PDF, Q g ð Þ may be derived if the velocity measurement resolves a sizeable portion of the dissipative range; i.e. PIV resolution well resolves the Taylor microscale. Given that the PDF is derived from the histogram of the occurrence of g, and that the measured variable in the above inequality is du g , it is straight forward to propagate the percent error in the instantaneous velocity of the PIV measurements to the uncertainty in the PDF. Our uncertainty of 2% in velocity translates to an uncertainty of 2.8% in du g . Given that the inequality is the only qualifying criteria, the uncertainty in the PDF may be achieved by perturbing the upper and lower limits of the inequality. To be conservative, we recalculated the PDFs by incorporating a 10% variation in the limits and found that the resulting PDFs with this additional 10% uncertainty in du g insignificantly influenced the shape of the PDF.   The PDF P(t p ) was determined by directly calculating the occurrences of t p for every occurrence of g using equation 2. Nondimensionalization parameters were calculated for each interrogation point and are listed in Tables 1 and 2.
Results: Dimensional Pdfs of g and t p for Flow Near a Bileaflet Mechanical Heart Valve Figure 4 presents the probability density function Q g ð Þ at points of interest (defined in Figure 3) during specific phases of the cardiac cycle. Upstream of the valve, the peak of Q g ð Þ (representing the mode of the fluctuating g) is consistently observed to be between 100-200 mm during forward flow. The smallest eddies captured in the distribution function are in the range 50-70 mm. The same characteristics are observed upstream of the valve during closure phase. However, during mid-diastole, there is a significant change in Q g ð Þ with respect to the varying location of points. In particular, except for the furthest location from the centerline, the locations closer to the centerline correspond to a significant leftward shift of Q g ð Þ. For these locations, the mode of g is in the range 80-90 mm, with the smallest g about 40 mm. Downstream of the valve, there is significant variation in Q g ð Þ characteristics as a function of the lateral position relative to the centerline during acceleration, peak, and deceleration phases ( Figure 4). The smallest g is between 40-50 mm and the mode of g between 90-200 mm. Positions located within the side orifice jets above and below the centerline do not show significantly different Q g ð Þ characteristics compared to the upstream forward flow characteristics.
The above shifts in Q g ð Þ characteristics clearly point to the slight leftward shift whenever increased turbulence is expected. For instance, locations slightly off from the centerline upstream of the valve will experience the high shear regions of the regurgitation jet during mid-diastole. Similarly, the leftward shifting of Q g ð Þ, also coincides with the shear layer locations downstream of the two leaflets consistent with the results in a more classical flow problem [36]. For all the locations, it appears that the smallest eddies are in the neighborhood of about 40 mm, which is consistent with the literature [12]. On the contrary, these small eddies appear to be relatively rare events with most of the dissipative eddies in the turbulent zones around 80 mm. This may be lower in locations where significantly greater turbulence exists. Figure 5 presents the normalized probability density function Q g=g O ð Þ for each of the positions of interest upstream and downstream at different phases of the cardiac cycle. g O is defined as g O~L Re {0:72 L where L is the local integral length scale, and Re L~k 1=2 L=n as discussed in Refs. [18,26] is a scale very close to . Also presented is Q g=g O ð Þ for the case of homogenous isotropic turbulence (HIT) from highly resolved direct numerical simulations [26] for comparison. Figure 5 shows very good agreement of the experimentally derived PDFs with the HIT expectation. The scatter around the HIT expectation may be attributed to experimental uncertainty as well as weak dependence of Q g=g O ð Þ on local mean shearing [20]. Specifically, we have shown that this weak dependence occurs when the Corrsin length scale, g c~S ET=S 3 À Á 1=2 , where S is the mean shear magnitude, approaches the mean shear-viscosity length scale, n=S ð Þ 1=2 [20]. Nevertheless, from a practical standpoint these results confirm the largely valid assumption of universal small scale structures despite the highly pulsatile nature of the turbulent flow past a complex device. This is particularly significant given the classical assumptions behind universality of turbulence are highly restrictive (i.e. very high Reynolds number, fully developed, equilibrium, stationary etc.). Figure 6 presents the raw PDFs of t p . The peak (mode) of the distribution for t p ranges between 5 and 20 dynes/cm 2 with the right tail extending to about 60 dynes/cm 2 . Tables 3 and 4 list the maximum observed t p corresponding to the minimum g at each position for all recorded phases for the cardiac cycle. To assess the universality of P(t p ), we examine the normalized probability density function P(t p =t p K ) (see Figure 7) for the different points of The normalized PDFs show strong data collapse indicating a practical universality which is expected based on the collapse observed for Q g=g O ð Þ. Note that the maximum observable t p was roughly 2t p K regardless of location of measurement or cardiac phase for the specific problem. The tale of P(t p =t p K ) drops exponentially around this magnitude. A theoretical limit to the maximum shear stress t p max may be estimated utilizing the smallest possible sub-Kolmogorov scale g S~L Re L {1 [26,28] to

Resolving Historically Inconsistent Hemolysis in Laminar and Turbulent Pipe Flow
The argument introduced to relate turbulence properties to the shear stress environment is a simple balance of energy dissipation. Thus, regardless of whether a flow is laminar or turbulent, isotropic or not, energy dissipated must highly correlate to the fate of suspended cells. In fact, energy dissipation was classically recognized as the hemodynamic parameter that dictates hemolysis based on simple physical arguments [38,39]. We test if this can resolve one of the paradoxical results published in Ref. [13] where hemolysis was measured in fully developed laminar and turbulent pipe flows while maintaining the same wall shear stress. Given that the total shear stress of the mean flows for both these cases are identical, the conventional thinking expects the same levels of hemolysis. However, turbulent cases produced significantly higher hemolysis [13]. To resolve the lack of an appropriate physical interpretation, we present a revised graph (Figure 8) that plots the measured plasma free hemoglobin (a measure hemolysis) as a function of total energy dissipation rate (in watts) in the pipe based on available information of wall shear stress and Reynolds numbers. Briefly, based on the wall shear stress magnitudes, length of the pipe, and its diameter the pressure drop can be calculated using a simple force balance. We also calculated flow rates from the given Reynolds number, viscosity, and pipe diameter information. The pressure drop and flow rate were multiplied to yield the total energy dissipation rate in watts. Figure 8 confirms that indeed, regardless of whether the flow is laminar or turbulent, the total energy dissipated produces a single monotonic functional dependence for hemolysis. The functional form with the best correlation was an exponential fit with R 2 = 0.91 compared to linear, and power fit models. Nonetheless, this result supports our fundamental assumptions and confirms that total energy dissipation is the fundamental parameter that dictates the mechanical environment of suspended cells.
One must note however that in the above pipe flow data, exposure time was not a relevant quantity as in all the experiments of ref. [13], the total exposure time was fixed. Thus the resulting trend shown in Figure 8 reflects the exponential growth in hemolysis due to the complex instantaneous energy dissipation ''histories'' of the turbulent cases. Further studies are necessary to develop appropriate blood damage models that take into account instantaneous energy dissipation rate histories and exposure times for pulsatile applications.

In the Context of a Cellular Environment
Cells are highly responsive to changes in their energetic environment, as demonstrated by the effectiveness of in vivo laser injury models [40]. Thermal energy is capable of inducing alterations to RBC morphology and can result in hemolysis in extreme cases [41,42]. Exacerbating the problem, a cascade of events can occur upon RBC lysis including platelet activation stimulated from adenosine diphosphate [43]. More direct studies of mechanical stress have demonstrated numerous effects on suspended blood cells, which can be combined with an energy balance equation, as presented here to determine the degree of cell activity. High shear stress is known to cause hemolysis and can further result in shear-induced platelet activation and shearinduced platelet aggregation, which has become known as SIPA, a feature that occurs independent of activation agonists. This process depends on the magnitude and the duration of shear stress [33]. More complex gradients in shear stress have been shown to promote activation-independent platelet aggregation involving the binding protein, von Willebrand factor (VWF) [44], a feature that may be important for the intermittent nature of turbulence. Platelets combined with VWF can result in aggregation while in suspension, independent of an adhesive surface substrate [45]. For sufficient shear stress, VWF can change from a globular confirmation to an extended chain [46,47], which can result in self-association, creating a structure that is very adherent to platelets and can create a physical barrier for other cells [48]. Shear stress in bulk flow or at a surface will also influence the inflammatory response. Activated platelets, alone, release a number of prothrombotic and proinflammatory molecules through a-granules and microparticles [49][50][51]. These processes may interact with mechanotransductive mechanisms in monocytes, which are highly responsive to membrane tension through cellular activation, signaling, and migration [52,53].
The current model may better predict the local shear stress environment for these processes and for dictating the function of binding proteins, platelets, RBCs, and monocytes, especially for the complex flow environment that can exist in medical devices and atherosclerosis [54]. As it becomes clear how cells respond to their mechanical environment, it will be increasingly important to develop unifying theories, such as the one presented in this paper to better determine the mechanical conditions for cells in a physiological or pathophysiological environment. It will further be critical for medical device development to consider the relationship between the cellular response and mechanics. For example, it has already been shown that self-association of VWF can result in acquired von Willebrand Disease in devices, including ventricular assist devices [55,56], a blood disorder that can be amplified by shear-induced receptor shedding in platelets [57].
In addition to predicting cellular functional responses to the local mechanical environment, it is also important to consider the influence on transport in the fluid flow. The complex biorheology of blood is highly dependent on the local shear environment [58]. The local environment influences the number of collisions between cells and the rotational nature of the cells, both of which are largely dependent on the hematocrit. Transport in the flow may be critical for thrombus growth rates and the movement of microparticles or cellular agonists [59,60].

Limitations
The proposed theory assumes that the distribution of eddies are not largely influenced by factors, such as the mean shear in the flow environment, an assumption shown to be reasonably valid for the mean shear at the points examined. However, there may be much stronger shearing in medical devices such as LVADs, where strong mean shear effects need to be accounted for to revise the distribution function of the dissipative scales of equation 3. In our probability density functions, we did see slight shifts, which may be caused by these influences, discussed in detail elsewhere [20,37]. However, the shifts were relatively minor in the current study.

Conclusion
We introduced a theory to predict the mechanical environment of cells in turbulent blood flows through physical arguments applicable to any in-homogenous turbulent flow. It is argued that for the case of physiological blood flow, the dissipative turbulent eddies are well in the Newtonian regime of blood rheology. We experimentally showed that the ''universal'' prediction of dissipative eddy distributions is valid even in a highly complex flow generated past a bileaflet mechanical heart valve, and that consequently the shear stress distributions experienced by blood cells must follow a practically universal distribution function. The universal form for shear stress acting on cells is presented and the theoretical maximum shear stress in a turbulent flow defined. Finally, the underlying assumption that local energy dissipation rate is the key to predicting the stress environment of cells is tested by demonstrating the high correlation between hemolysis and energy dissipation in pipe flow regardless of laminar or turbulent flow regime.