Partitioning of red blood cell aggregates in bifurcating microscale flows

Microvascular flows are often considered to be free of red blood cell aggregates, however, recent studies have demonstrated that aggregates are present throughout the microvasculature, affecting cell distribution and blood perfusion. This work reports on the spatial distribution of red blood cell aggregates in a T-shaped bifurcation on the scale of a large microvessel. Non-aggregating and aggregating human red blood cell suspensions were studied for a range of flow splits in the daughter branches of the bifurcation. Aggregate sizes were determined using image processing. The mean aggregate size was marginally increased in the daughter branches for a range of flow rates, mainly due to the lower shear conditions and the close cell and aggregate proximity therein. A counterintuitive decrease in the mean aggregate size was apparent in the lower flow rate branches. This was attributed to the existence of regions depleted by aggregates of certain sizes in the parent branch, and to the change in the exact flow split location in the T-junction with flow ratio. The findings of the present investigation may have significant implications for microvascular flows and may help explain why the effects of physiological RBC aggregation are not deleterious in terms of in vivo vascular resistance.


Results
The experimental set-up shown in Fig. 1 was utilized in the study. The sample was perfused through a microchannel with dimensions of W = 100 μm (width) and D = 40 μm (depth) by applying pressure to the sealed inlet reservoir. The distribution of the flow in the two daughter branches (flow split) was altered by changing the height difference between the two outlet reservoirs. The origin of the global coordinate system was placed at the intersection of the central axis of the parent and daughter branches and the coordinates were normalised with the channel width (x* = x/W and y* = y/W). Figure 2 shows representative raw and processed images for the non-aggregating (PBS) and aggregating (D2000) blood cases, for flow ratios of Q* = Q D /Q P = 0.5 and 0.2 respectively. The schematic in Fig. 3 illustrates the definition of the main quantities of interest in this work; the area of a detected structure in an image (A, in μ m 2 ), the normalized area of a detected structure (A*), the mean normalised area size in one image ( ⁎ A ), and the mean normalised area size from all images (ˆ⁎ A ). The normalising parameter is the characteristic area of one RBC (A c , in μ m 2 ). The size of aggregated structures present in the flow (i.e. the detected areas) was estimated using an edge-detection image processing technique 24 . The technique takes advantage of the fact that isointensity regions develop in the image as a result of the RBC aggregation process. Further details on the experimental parameters and methodology are provided in the Methods section.
Aggregate size distribution and haematocrit. Figure 4a compares the normalised size (A*) distributions of the structures detected in the entire channel (number of images n = 400) for aggregating (D2000) and non-aggregating (PBS) cases respectively, and for a flow ratio of Q*~0.5. It should be noted that the distributions are left-truncated at 0.5 for the purpose of noise reduction. A* = 1 indicates the characteristic area of one RBC. However, red blood cells may adopt various configurations as they flow through the channel, altering the Regions of interest are shown as white, dashed rectangles. The velocity field and streamlines derived from an aggregating case with Q*~0.13 is shown in the right panel. The vector and streamline density is reduced by half for clarity of presentation. The red lines are the lines that separate the flow to the different daughter branches. The dimensionless quantities y* and y* are the x and y coordinates normalised by the channel width (100 μ m). projected, and hence detected, RBC edge area in the images and resulting in the A* distribution shown in Fig. 4(a) for the non-aggregating case. The detected structures are significantly larger in the aggregating (D2000) case exhibiting a mean A* value of 1.90[0.50′, 1.91] (k = 73164, n = 400). Mean values are reported with their lower and upper bounds of the 95% confidence intervals, with the primed values indicating truncated data (k is the number of structures). A similar size distribution with a mean A* value of 1.88[0.50′, 1.90] can also be observed at a lower flow ratio (Q*~0.1) for the aggregating case ( Fig. 4(b). These mean values suggest an almost twofold increase in the mean area size of aggregates compared to that of one RBC. It should be noted, however, that in the aggregating cases the area range of one RBC (minimum < A* < 3) does not necessarily imply structures of dispersed RBCs, since such areas could also be parts of aggregates. Excluding this area range from the distribution could provide information about larger aggregate characteristics; the mean area A* increases to ≈ 5.59[3′, 5.68] (k = 12378) in this case.
In order to examine the effect of local haematocrit on the distribution of A*, the spatial distribution of the detected structures across the parent and daughter branches is plotted in Fig. 4(c) and (d) respectively for the PBS case. These distributions show the location of all the detected structures in the channel irrespective of size. Corresponding haematocrit profiles estimated from the image intensity, H(I*), are also superimposed on the same figures. These intensity-dependent haematocrit profiles were developed in Sherwood et al. 20,22 and represent linear (H l , dashed lines) and non-linear (H n−l , solid lines) relationships. In the parent branch (Fig. 4c) the distribution of detected structures generally matches the estimated haematocrit profiles. However, this is not the case in the daughter branches. This can be attributed partly to the assumed linear relationship between haematocrit and image intensity in the H l (I*) profiles, and to a small sensitivity of the image processing technique to the cell concentration (local haematocrit) 24 . Additionally, the calibration of the non-linear function H n _ l (I*) was performed in a microchannel of slightly greater depth (50 μ m cf 40 μ m in the present study). It should be noted that no significant dependency of the spatial distributions on structure size was noted; i.e. spatial distributions plotted for different structure size ranges (for instance 1 < A* < 2, or 2 < A* < 3) remained qualitatively similar.
Effect of flow ratio. Our previous work has shown that in the T-junction, the RBC concentration (haematocrit) decreases with the decrease of Q* 20 . For the non-aggregating case no dependency of the detected edge areas on flow ratio Q*, and therefore haematocrit, should be observed. This is illustrated in Fig. 5(a), by plotting the ensemble average ˆ⁎ A (n = 400) against Q*. The pooled standard deviation σˆ⁎ pA (equation 3 in Methods section) is shown in Fig. 5(a), and the standard deviation σˆ⁎ A is shown in Fig. 5(c). The ˆ⁎ A values in Fig. 5(c) were normalised by the values calculated in the parent branch and noted as ˆ⁎ A n . The magnitude of the pooled standard deviation seen in Fig. 5(a) relates to the fact that the random orientation the dispersed RBCs assume within the channel results in a range of detected structure areas. The low standard deviation values observed in Fig. 5c indicate the consistency of the processing and the reproducibility of the technique.
The flow ratio Q*, however, has a distinct effect on the distribution of aggregated structures in the channel. This becomes apparent when the relationship ˆ⁎ A -Q* ( Fig. 5(b) and (d)) is compared to that of the PBS cases.

. ˆ⁎
A values tend to 1 at higher Q* ratios as expected. The preference of the larger aggregates for higher flow ratio branches is depicted more clearly in Fig. 5(e), where a size-flow parameter, F* A* = Q*A*, is plotted against Q*; above approximately Q* = 0.3 the size-flow parameter F* A* deviates from the F* A* = Q* line, indicating an increase of aggregate size relative to the mean aggregate size in the parent branch.
The decrease of ˆ⁎ A n , however, below Q* = 0.2 is counterintuitive considering the very low shear conditions in the daughter branches for these flow ratios, which should promote RBC aggregation. The cause of this behaviour can be traced to the structural characteristics of blood and the flow characteristics in the channel. As already mentioned, one of the main characteristics of the flow in the T-junction region is the flow split which is determined by the incoming flow (from the parent branch) and the pressure difference in the daughter branches. The flow split location in the parent branch, i.e. the boundary separating the flow entering each of the two daughter branches, is determined by the separating streamlines in the T-junction region (please see Fig. 1). The flow split location (denoted as ⁎ x fs ) is expected to depend on the flow ratio Q* and hence to deviate from the channel centreline for asymmetric flow splits. ⁎ x fs indicates the location of the flow split line on the x* axis. Partitioning of aggregates in the channel. In order to examine the preferential location of the detected aggregates in the channel, spatial distributions for two distinct A* ranges (3 < A* < 5 and 7 < A* < A* max ) were produced and shown in Fig. 6(a) for the parent branch and for all flow ratio cases examined in this study (for Q~0.35 ± 0.5 μ l/h). The particular ranges of A* were chosen for illustrative purposes. Also, in order to exclude the sizes corresponding to one RBC, sizes up to A*~3 are excluded. The distributions in Fig. 6(a) illustrate that the larger structures (A*~7 to maximum size) tend to concentrate in the centre of the channel due to lower shearing forces therein, resulting in an aggregate depleted region near the wall. In addition, when larger aggregates form there is less cell-cell interaction and therefore the radial dispersion of the structures is decreased 30 ; the wall-induced lift forces have a stronger effect in concentrating the aggregates towards the flow centreline. One implication of the concentration of the large aggregates around the flow centreline is a preference of these structures to flow towards the branch with the highest flow rate. This is apparent in Fig. 6(b) which shows the distribution of three different sizes of aggregates (A*~3 ± 0.02, 5 ± 0.02 and 7 ± 0.02) in the channel, for Q*~0.13. The flow separating streamlines are also shown in the figure to aid the discussion in the following sections.
Effect of flow ratio on aggregate partitioning. To illustrate the preferential concentration of larger aggregates towards the flow centreline in the parent branch and the resulting distinct aggregate depleted region near the wall indicated in Fig. 6, the spatial distribution of an arbitrary selected aggregate size (A*~4) in the parent branch, is plotted in Fig. 7(a). The arrow indicates the location of the mean value of A*, calculated from the data falling below the 5 th percentile of the distribution (right side); the location on the x* axis is denoted as ⁎ ⁎ x A . In the region between ⁎ ⁎ x A and the parent branch wall aggregates of size A* around 4 are generally absent. This region, together with the separating line for a given Q*, influences the partitioning of aggregates in the bifurcation and hence their distribution in the daughter branches. It is expected that for | | > | | ⁎ ⁎ x x fs , aggregates will flow into the adjacent daughter branch.
The flow-split location ⁎ x fs for all Q* values is shown in Fig. 7(b). The x* A* locations in the parent branch, are plotted in the same graph as lines parallel to the x-axis, since they are independent of Q*. These are estimated  from 400 images for each of the 7 distinct aggregate sizes (i.e. A* = 1 to 7) shown. The sizes of aggregates that will flow in a branch having a specific Q* value can be estimated from Fig. 7(b): for a specific Q*, the A* sizes found in the shaded area defined by the estimated low split locations are excluded from the particular branch, and only the A* sizes above a specific x* f-s value (black-dots) will be transported to the relevant branch. For example, at Q*~0.2, aggregate sizes corresponding to values above the flow split location x* fs for the given flow ratio, i.e. A* ranging from 1 to 5 are found. Aggregates in that size range will be transported into the branch with Q* = 0.2, whereas aggregates of all sizes will be transported to the opposite branch with Q*~0.8. Note that at Q*~0.5 the flow split line coincides with the centreline of the channel flow including all aggregate sizes.

Discussion
The analysis of the non-aggregating cases (PBS samples) served primarily as a technique validation exercise. The main point of interest in these cases is the spatial concentration of structures in the channel, which has been studied previously 19,20,22 . Other interesting features of this flow, such as the orientation of RBCs and/or their frequency of rotation in the flow, could not be resolved at the haematocrit employed in the study. To study the details of the RBC motion at physiological haematocrit, either special treatment of the samples is needed 30 , or computational models can be employed 31 , although such models require rigorous validation with experimental data.
The analysis of the aggregating samples (D2000 data) revealed valuable new information regarding the partitioning of the aggregates in the T-shaped bifurcation and their distribution further downstream. In Fig. 5(d) it is observed that ˆ⁎ A n increases above 1 (i.e. exceeds the parent branch ˆ⁎ A value) for Q* above ~0.2. This can be attributed to the reduction of the shear forces as a result of the lower flow rates in the daughter branches compared to the parent branch. Furthermore, after partitioning at the bifurcation aggregates are forced to one side of the channel causing a higher compactness and rate of interaction between them. Although the exact shear rate distribution in the channel cannot be obtained without knowledge of the flow field in the depth direction, it is reasonable to expect the mean shear in the daughter branches to scale with the mean velocity and hence with Q* (as the parent branch flow rate was approximately constant for all cases). An approximate measure of the level of shear is given by the ratio of the mean velocity (U), to the depth of channel, i.e. γ  = U μ m s −1 /40 μ m. The magnitude of this 'pseudo-shear' in the parent branch is approximately 8 s −1 , which is considered moderate to low in terms of its influence on aggregate dispersion. In the daughter branches the pseudo-shear scales with Q* by definition, reaching very small values for the lower flow ratios.
The first novel piece of information extracted from Fig. 6(a) is that as the size of aggregates increases so does their distance from the side-walls of the channel. This behaviour is illustrated in Fig. 7b, where the distance x* A* was plotted for various aggregate sizes (A* = 1 to 7). The development of aggregate structures of various sizes is mainly attributed to the flow conditions, and more specifically the local shear rates, which determine the magnitude of the shear forces in the channel. The dependency of the RBC aggregation phenomenon on shear rate is well established in the literature [32][33][34] . However, estimating the aggregate size distribution at physiological haematocrits is not a trivial task; hence it is not surprising that little such information exists in the literature. To the best of our knowledge no detailed analysis of aggregate characteristics in geometries more complex than a straight microchannel and at normal haematocrits has been conducted hitherto.
The behaviour of x* A* in Fig. 7(b) implies that, in addition to the haematocrit variation in the cross-flow direction (illustrated in Fig. 4(c) and (d) and analysed in Sherwood et al. 19,20,22 ), the size of the aggregates varies also. This implies that there are regions adjacent to the side-walls of the channel that are depleted of aggregates of certain sizes. Indeed, in Fig. 6(a) a clear region depleted of aggregates can be observed, although to term this an aggregate depleted layer would be an oversimplification, as the region is aggregate size-dependent. The implications of the existence of these regions can be seen in Figs 5(d) and 6(b), where a preferential distribution of aggregates is evident in branches with low flow ratio. Although the dependence of blood viscosity on RBC aggregation, in addition to other factors such as haematocrit, plasma viscosity, etc. is well established in simple geometries, the additional structural characteristics of aggregating blood found in this study provide significant information towards a better estimation of local viscosity profiles in more complex environments.
In a previous study investigating the local viscosity characteristics based on the spatial variation of RBC concentration we 20 utilised the haematocrit-dependent viscosity model of Pries et al. 21 to create viscosity profiles in the microchannel daughter branches for aggregating and non-aggregating cases, and illustrated the significant effect of RBC aggregation in this aspect. Local aggregation information, not included in the models of the aforementioned studies, was incorporated in another work by Kaliviotis et al. 35 , for aggregating blood flows within a parallel plate geometry. It was shown that significant spatial variation in the viscosity exists due to aggregate formation in the flow; the viscosity model used in that work was developed to account for various aspects of RBC aggregation including aggregate network formation 36,37 . The presence of the formed aggregates in the flow and the estimated distributions are expected to further affect the local viscosity of the fluid, since the mechanical characteristics of such structures are considerably different to those of RBCs in the dispersed state, even at similar local haematocrits. This is well established in the literature through studies such as the seminal work of Chien et al. 38,39 .
Accurate estimation of the local viscosity distribution in an aggregating sample, however, requires characterisation of the shearing conditions in the channel, which in the present analysis is restricted due to the lack of three dimensional velocity data, the aspect ratio of the channel and the depth of field of the microscopy system 24 . The RBC aggregates, however, could provide information on the extent of shearing they experience in the flow, as they are very sensitive to the local shear forces; in that sense, the aggregate size distributions presented in Fig. 4(a) and (b) imply that a range of shearing rates exists in the flow.
The aggregate size partitioning behaviour seen in Figs 5,6 and 7 may help explain why the effects of RBC aggregation are not as deleterious as expected 15 in terms of in vivo vascular resistance; the smaller aggregates entering the low flow rate (and therefore low mean shear rate) branches, in combination with the lower haematocrit in these regions 20 , imply that the bulk viscosity in the branch will be kept low as well. This is a significant implication as it is known that blood exhibits yield stress 40,41 at very low shear rates, which would have a negative impact if occurred in the low flow branches. Furthermore, since oxygen perfusion takes place mainly in the capillary level, the large CDL, and the low concentration of the aggregated cells in larger scale branches with low flow rates, as found in the present study, would serve the purpose of flow promotion. It should be noted here that the concentration of cells and aggregates in the low flow branch is strongly skewed towards the wall, i.e. in the lower velocity, but higher shear, regions of the flow. This may contribute to reducing the size of aggregates in the branch.
The present study considers a T-junction geometry (i.e. 90 degree symmetric bifurcations) which typically exists in the microcirculation 42 . To extent the findings of this study to other bifurcating geometries, such asymmetric T-junctions etc., the two main aspects that were found to affect the partitioning of aggregates should be considered: the flow split location x* fs and the aggregate depleted region, defined by x* A* in the parent branch. The aggregate depleted region is independent of the flow ratio and should depend on the aggregation intensity of blood. The flow split location, as shown in studies using bifurcations other than the T-junction 43 , seems to depend on the flow-split ratio in the outlets (Q*) only. This implies that the current findings could be extrapolated to geometries other than the T-junction. However, further experimental data might be required to verify this.
Another important flow parameter in the microcirculation is the haematocrit, as it largely affects the bulk and local viscosity of blood 20,21 . The haematocrit used in the present study (25%) is considered physiological for the length-scales studied 44,45 . However, the haematocrit may be altered in various pathologies. A change in the haematocrit -while flow conditions and RBC aggregation intensities remain the same -can affect the shear stress distribution in the flow, and subsequently the extent of aggregation, as the later depends on the shearing forces developed in the flow. The present study has shed light into important and previously unexplored structural phenomena of aggregating blood in the microcirculation. Nevertheless, further work is required to fully address both the influence of the bifurcation type and asymmetry, and the influence of haematocrit on local aggregation characteristics.

Conclusions
RBC aggregation affects many aspects of blood flow at low shear conditions, including CDL, local haematocrit and local viscosity distributions, which have been demonstrated in previous studies by the authors 19,20,22,24 using rectangular microchannels. In the present work additional structural characteristics, i.e. the size distribution of the formed aggregates in a T-shaped bifurcating flow, were quantified.
The significance of characterising aggregation at a local level is illustrated by the novel findings of this study; the analysis of the aggregate transport to the different daughter branches revealed a counter-intuitive behaviour between aggregate size and flow conditions, i.e. small aggregates found in branches of very low shear conditions in contrary to the current understanding on the aggregation/shear relationship. This phenomenon was attributed to the non-uniform distributions of aggregates in the parent branch, which in combination with the existence of a near wall region depleted of aggregates of certain sizes, affect the partitioning of aggregates for any given flow ratio, and hence the aggregate size distribution in the daughter branches. This new information on aggregate structure size and distribution could improve the estimation of local blood viscosity in the microscale and enhance further our understanding of microvascular flows. It may also help explain the contrast between in vivo and in vitro observations regarding the effects of RBC aggregation on vascular resistance or viscosity.

Methods
The details of the experimental apparatus, sample preparation and flow measurement methodology can be found in Sherwood et al. 19,20 whereas that of aggregate size characterization in Kaliviotis et al. 24 Only a brief description is provided below. Sample preparation. Sample acquisition and all procedures involving human participation were performed with the approval, and in accordance with the relevant guidelines and regulations set by the Southeast London Ethics Committee (ref: 10/H0804/21). Blood was collected from healthy volunteers into vacuum tubes (BD) preloaded with 1.8 mg/ml EDTA after obtaining informed consent. The RBCs were washed twice in Phosphate Buffered Saline (PBS), centrifuged at 3000 rpm, and suspended in PBS containing D2000 (5 g/l). The haematocrit was adjusted to 25% by volume as appropriate for the microchannel dimensions used 44,45 . For consistency all experiments were conducted with a single sample.
The flow system. A microchannel, with a width of W = 100 μm and depth D = 40 μm, was fabricated from SU8 using photolithography (Epigem, Redcar, UK), and is shown in Fig. 1. The sample was perfused through the microchannel by applying a pressure to the sealed inlet reservoir, relative to the open outlet reservoirs. The pressure in the inlet reservoir was controlled with an in house pressure control system using an actuated needle valve and a compressed nitrogen source. In order to minimise effects of RBC sedimentation, which is particularly pronounced for aggregating samples, the fluid in the inlet reservoir was continuously mixed with a magnetic stir bar, except when acquiring data, and all tube lengths were kept to a minimum.
For the cases studied, the inlet reservoir pressure was set to a single value, resulting at a mean velocity of 320 μ m s −1 in the parent branch. The distribution of flow between the two daughter branches (flow split) was altered by means of hydrostatic pressure difference. The pressure difference was achieved by independently adjusting the height of the outlet reservoirs using micrometer stages. The flow split for a given acquisition case was then calculated using the PIV data, as described in the following paragraph. Between acquisitions, the channel was perfused at a high flow rate in order to ensure uniform hematocrit throughout the channel and system. Following, the pressure was reduced to the desired level and 20 s were allowed for aggregation to reach a steady state before acquisition proceeded 20 .
Scientific RepoRts | 7:44563 | DOI: 10.1038/srep44563 The origin of the global coordinate system is shown in Fig. 2a at the intersection of the parent and daughter branches. The coordinates are normalised with the channel width, x* = x/W and y* = y/W, respectively, and thus the parent branch width spans from x* = − 0.5 to 0.5 and the daughter branches from y* = − 0.5 to 0.5.
One T-junction channel was used throughout the experiment and 2000 images were acquired for each flow ratio Q*. The same blood sample was used for consistency at each Q* (in the range 0-1). Microscopy and microPIV system. The flow system was mounted on the stage of an inverted microscope (Leica DM ILM, Germany), with the focal plane set to the centre of the channel. A halogen light-source was used for illumination and images were acquired using an IDT X3 CMOS camera (Tallahassee, USA) at a frequency of 125 Hz. After pre-processing the images for alignment with the particle image velocimetry (PIV) interrogation windows, a final image size of 1216 × 700 pixels at a spatial resolution of 0.65 μm/pixel was obtained. Multi-pass ensemble averaged PIV processing was carried out on each of the data sets, providing a final window size of 8 × 8 pixels and a vector spacing of 4 pixels (2.6 μ m) using JPIV software (www.jpiv.vennemann-online.de/). A normalised median test was utilised to identify the invalid vectors, which were replaced with the median of the surrounding vectors. For the estimation of the flow rate for each branch, spatially averaged velocity profiles were acquired in the regions indicated by the region of interest in Fig. 1. The flow ratio was defined as the ratio between the flow rate in the daughter and the parent branch, Q* = Q D /Q P, and was estimated from the average velocities (U) measured with PIV, according to Q* = U D /U P , as the cross-sectional area was the same in all branches.
The velocity fields were produced using the JPIV software. The velocity vectors are mean values from 2000 processed images. The data were further processed in Matlab to derive flow streamlines and in particular the separating streamlines that provide information on the flow partitioning at the bifurcation as a function of flow ratio (see Fig. 1). The velocity field and streamlines in Fig. 1 were derived from an aggregating case with Q*~0.13. The vector and streamline density is reduced for clarity of presentation. The red lines are the lines that separate the flow to the different daughter branches, and the point where these lines originate from in the parent branch is defined as the flow-split location in the x* direction (noted as x* fs ). The flow split location was derived for all Q* cases analysed.

Image based detection of aggregates.
The method for detecting RBC aggregates in the acquired images was developed by taking advantage of the fact that connected RBCs (often termed rouleaux) result in continuous iso-intensity patterns, which can be detected with edge detection algorithms. In Kaliviotis et al. 24 it was shown that these continuous areas of similar intensity increase almost linearly with the number of RBCs per aggregate, with relatively little deviation, for small aggregates of up to 15 RBCs per aggregate.
The detection algorithm is detailed in Kaliviotis et al. 24 and briefly described below. It comprises the main image processing stages: • Image pre-processing to improve image quality; these included: intensity transformations for improving the global contrast of the image 46 and a correction of uneven illumination using a Gaussian function. 22 • Local intensity gradients (∇I i ) identification: the magnitude of the local intensity gradients was detected using an algorithm based on the Prewitt method 46 . Four intensity differences per distance (Δ I/Δ S) were computed from all eight opposite neighbouring pixels in a 3 × 3 pixel convolution window and then the magnitude of the local intensity gradient was calculated by applying spatial convolution: The magnitude of the local gradient (∇ I i ) was normalised by the mean intensity of the convolved region (I con ): −1 ]. This was necessary because the gradients produced by aggregates are influenced by the background intensity; this dependency is complex to analyse when considering the whole gradient spectrum, however, when focussing only on the edges of RBCs it tends to be non-linear following a parabolic behaviour 24 .
• Blurring correction: small areas at the extreme ends of the daughter branches, were slightly blurred. Although this blurring does not have a significant effect on the cross-correlation process of the micro-PIV technique (as the PIV measurement is depth saturated), it does affect the local gradients ∇ ⁎ I i within those areas. Therefore a correction procedure was implemented for the blurred regions, which was based on creating a mean con- producing an image with pixel values equal to 1, for the unaffected regions, and pixel values lower than 1 for the blurred regions (~0.87 was the minimum value). The affected gradients in the blurred regions were corrected using the following expression: • Edge detection: The points in which the local gradients were above a certain threshold value were defined as edge points and converted to white (assigned a value of 1). Pixels with values below the threshold were assigned a value of 0. Calibration of the gradient threshold value was performed in the non-aggregating cases (PBS) and used blindly in the rest of the aggregation cases. It involved the identification of a threshold value that would result in a mean edge area (interconnected white pixels) equal to a characteristic RBC area A c . This characteristic RBC area was defined based on the fact that different configurations of one RBC could be seen in the images due to the rotation and deformation of the cells in the flow; the estimated edge size of one RBC when viewed from the flat side, i.e. its maximum edge size, was estimated at ~45 μ m 2 , whereas the edge area of one RBC when seen laterally (the minimum value) was estimated to be ~10 μ m 2 . A conservative value of A c = 15 μ m 2 was set based on the fact that the various limitations imposed by the nature of the imaged data (dense suspension, overlapping of cells, relatively low resolution, etc.), affected the overall quality of the images and prevented, to a large extent, complete RBC edges (e.g. in their flat position) to be resolved. • Noise elimination: The detected edges (defined by the interconnected white pixels) that were smaller than 20% of the mean RBC edge area were discarded as noise. • Statistical analysis: Statistical information about the identified edges that appeared in the processed image was extracted for further analysis; this included number, area and location of aggregates.
The methodology has been validated against other image processing and electrorheology techniques 47 in Kaliviotis et al. 24 and its limitations have been discussed therein 24 . Figure 2 shows representative cases of processed images; the original image and the detected edge areas of a PBS case for Q*~0.50 is illustrated in 2(a) and 2(b), whereas the result for the aggregating case (D2000 sample) for a Q* ~0.2 is shown in Fig. 2(c) and (d). In 2(d) small areas (50% of the mean RBC edge area, ~8 μ m 2 ), have been omitted for clarity of presentation of the larger structures. The size of the structures in Fig. 2(b) and (d) is indicated by using a colour-scale for a qualitative comparison; the sizes have been normalised relative to the largest size in the D2000 case. As it can be seen in Fig. 2(b) there is a narrower size distribution (mostly dark coloured structures corresponding to smaller sizes) compared to the D2000 case in Fig. 2(d) (lighter colours corresponding to larger structures).
Statistical analysis. The measured and calculated quantities (pixel intensity, velocity magnitude, edge area, etc.) are typically expressed in terms of spatially and/or temporally averaged values. The normality of the data distributions was assessed with the one-sample Anderson-Darling test, using Matlab (Mathworks) software.
For normal distributions, the mean values (MV) and standard deviation (SD) are presented as MV ± SD. For exponentially distributed data, which are observed for the spatially varying edge areas in one image, the mean values (MV) and their lower (LB) and upper bounds (UB) of the 95% confidence intervals were obtained using a least squares fitting function 48  The detected edge area of a single structure in an instantaneous image is denoted by A and when normalised by the characteristic area of one RBC (A/A c ) it is denoted by A* (see schematic in Fig. 3). The mean value of k detected structures in one image, normalised by the characteristic area of one RBC (A c ), is denoted by ⁎ A with standard deviation σ ⁎ A . Note that for the exponential distributions of the spatially-varying A*, which is the case in the present work, σ ⁎ A is equal to ⁎ A by definition. When the time-varying mean value of ⁎ A (denoted by ˆ⁎ A ) is calculated from a number (n) of images two relevant quantities regarding the dispersion of the data are relevant: a) the standard deviation of ˆ⁎ A , which is derived from the time-varying data and is denoted as σˆ⁎ A (note that σˆ⁎ A is calculated from normally distributed data), and b) the pooled (combined) standard deviation of ˆ⁎ A , which is derived from the spatially-varying data in individual images and is denoted as σˆ⁎ pA ; the latter is derived from exponentially distributed data, as explained above and provides a measure of the average dispersion of the data in the images. A sample of 400 images (n = 400) was analysed (equispaced out of the 2000 images acquired) for all cases and the number of structures (k) is indicated in each different case. The abovementioned quantities are calculated from the variances:  The constant α in H l determines the actual haematocrit, however, in the present case it was adjusted for qualitative comparison purposes with the spatial distribution of A*; this was necessary due to the fact that the number of detected structures was not scaled for haematocrit. The values of the parameters a and b were defined in Sherwood et al. 20,22 after calibration (a = 0.685, b= 9.244) in a 50 × 50 μ m rectangular microchannel.
Limitations and error. The limitations of the image processing technique have been discussed in detail in Kaliviotis et al. 24 The main factor affecting the image processing is the depth of focus which covers the entire depth of the channel, and therefore, overlapped cells appear in the image. This issue was addressed by calibrating the parameters for conservative processing, that is utilising the higher range of the image intensity gradients for the definition of edges and structures. This results in a relatively low density of detected structures, compared to the actual cells in the flow as already mentioned, and an underestimation of the aggregate sizes in the flow. Hence the main source of error in the data presented in this study is the parameter A*. Based on the validation process described in Kaliviotis et al. 24 the error in the estimation of structure size A* might be of the order of 15% (estimated from the correlation curve of A* to actual aggregate size measured by inspection) with considerable deviation for aggregates consisting of more than 15 RBCs. Assuming that this error is systematic in the analysis it can be concluded that the trends shown in the present work are representative of the actual phenomena taking place in the flow configuration. Note that the standard deviation of the data in Fig. 5(d) is of similar order to the expected error.