Exploring the potential of blood flow network data

To gain a better understanding of the role of haemodynamic forces during the development of the cardiovascular system, a series of studies have been reported recently that describe flow fields in the vasculature of model systems. Such data sets, in particular those reporting networks at multiple stages, mark a transition in the focus from single blood vessels to large parts of vascular networks. It becomes possible to investigate the behaviour of a blood vessel in the context of its surroundings, rather than as an isolated entity. In this study, a framework is presented that facilitates the analysis of such data sets. The blood vessel data is represented as a graph, with each node connected by a vessel segment with known properties. Using this framework the pressure distribution and other parameters of interest can then be estimated. Two examples are given that make use of this scheme: (1) a method to detect and reduce measurement errors in the network and (2) a method that allows the testing of various haemorheological models. For both examples a proof-of-principle result is shown.


Introduction
The role of fluid mechanics in the development of vascular networks has received considerable attention in recent years.Haemodynamic forces, such as the local wall shear stress, have been suggested as being essential epigenetic factors during cardiovascular development [7,14,19].This has led to the vascular remodelling paradigm, as shown schematically in Fig. 1.The core idea here is that the flow through a network leads to continuous changes in its structure, in order to maximise a particular network function (e.g.nutrient transport), while at the same time minimising energy expenditure (e.g.due to viscous losses).This continuous remodelling by flow is thought to be able to convert an initially unstructured network (or 'plexus') into an efficient branching network.This behaviour has been observed, for instance, in mouse embryos [16] and chicken embryos [13,14], both common model systems to study human cardiovascular development.The advantage of these model systems is that they are also accessible for mechanical or chemical intervention, so that deviations from normal development can be observed (e.g. to see the effects of ligating a vessel).
While the requirement for flow for the proper development of the cardiovascular system is undisputed, the precise role is still poorly understood [19].The endothelial cells that line the vessel walls are well known to respond to haemodynamic forces [4], but how these local events can lead to global topological changes is still a mystery.This can be rephrased as: how can a complex, efficient vascular network emerge from simple, local rules?These local rules can ostensibly be based only on information that is available to the endothelial cells, e.g.wall shear stress, pressure or concentration of a certain component.
To better understand vascular remodelling, detailed experimental data sets are required.Until a few years ago, data was solely based on imaging, in order to obtain the network topology.Occasionally, the studies are supplemented with qualitative flow observations, for instance using ink visualisation.In recent years, however, rapid developments in non-invasive measurement techniques have made it possible to obtain simultaneous measurement of the shape and functioning of networks [10,12,25].These studies provide detailed flow information in a relatively large field-ofview.An example of one of the capabilities of these modern approaches is shown in Fig. 2, based on data obtained by Kloosterman et al. [13].In their approach, which is based on in vivo microscopic particle image velocimetry [20], a sequence of images is first obtained using digital cameras that document the motion of tracers (here erythrocytes).The local displacement is subsequently estimated by crosscorrelation of small regions of the total image sequence.The end result is a detailed velocity field, with a resolution of approximately 10 Â 10 lm 2 and a total field of view of typically 3 Â 5 mm 2 .Here only the time-averaged field is shown, but the pulsatile nature can also be retrieved [20].Note that implicitly, it is assumed that these velocity fields represent the center-plane of the more-or-less two-dimensional extra-embryonic vitelline network; volumetric measurements are also possible, by stacking different measurement planes [18].While these modern techniques provide an unprecedented level of detail, they also present a practical problem: how to analyse the vast amount of two-, three-or four-dimensional data that is available, especially if also many networks are obtained?A successful approach has been to make use of network models (in particular graphs, see Sect. 2), which greatly simplifies the handling and interpretation of the data.However, the majority of vascular graph studies are based solely on image data (and thus topology only).In the measurements such as shown in Fig. 2, also flow information is available.In this paper, the focus will be on how this additional information can provide various exciting opportunities to study haemodynamic phenomena.Two examples are given: (1) correction of flow data based on conservation of mass and (2) in vivo testing of haemorheological models.The chicken embryo data provided by Kloosterman et al. [13] will be used throughout this study, but the approach can readily be applied to other data sets.

Network approaches
Flow networks (or hydraulic circuits) occur in many engineering applications, such as water supply networks, blood vessel networks and microfluidic devices.Finding the complete, complex flow field in such large scale networks-e.g. using computational fluid dynamics-is generally impractical, but often also unnecessary.To make the analysis of these networks tractable, some important simplifications are often made: the network is represented as a onedimensional system, consisting solely of nodes and branches.Nodes only serve to distribute flow among the connected branches and impose no resistance to Fig. 1 The vascular remodelling paradigm: an initially unstructured network remodels into a hierarchical network due to a continuous process of remodelling, guided by haemodynamic forces the flow.The branches, on the other hand, represent the hydraulic resistances.These resistances are lumped parameters describing all head losses for the given branch geometry and flow conditions.For the network, two rules can be formulated, analogous to Kirchhoff's laws for electric circuits: (1) mass is conserved at each node, (2) the pressure drop along a closed loop is zero.The analogy between electric and hydraulic circuits also implies that similar analysis tools are available to find the flow distribution.The classic example is the method introduced by Hardy Cross, which allows one to efficiently find the flow distribution by hand, for particular pressure boundary conditions [3].
Over the years, there have been many studies that applied the network approach to the cardiovascular system.One of the pioneering works is the publication by Lipowsky et al. [15], who applied the network approach to data obtained in the cat mesentery.For the given network topology, they obtained the pressure and flow distribution in the network.Notably, they only reported statistics averaged over all segments of a given diameter; in other words, they did not consider the 'context' of a certain blood vessel.
A similar network approach has been used, among others, to study the vasculature of the retina [5], coronary blood flow of a pig [11], transport of oxygen in arteriolar networks in a hamster cheek pouch [21], and cerebral blood flow [23].Pries et al. [22] included microrheological effects, such as the Fa ˚hraeus-Lindqvist effect (i.e. to account for the variation of viscosity with vessel diameter).Van den Wijngaard et al. [26] extended the method to three dimensions to study the coronary circulation.Apart from these applications based on existing physiological data, the network approach has also been used in more theoretical studies, for instance to study the emergence of patterns in initially hexagonal blood vessel networks [6].Similarly, by including mass transport and vessel compliance, Boas et al. [2] were able to show dynamic behaviour in a relatively simple model.
Note that virtually all applications of network analysis, the strategy was the same: for a given network layout (either from imaging data or prescribed, e.g. as a hexagonal lattice) the authors find the pressure and flow rate distribution.From the latter, derived quantities (local velocity, wall shear stress) can be obtained.Pressure and flow are found simultaneously using the common assumptions and solution method (see next subsection).However, in many of the new experiments the flow data is actually available, see e.g.Fig. 2.This means that less assumptions have to be made when finding the pressure.Furthermore, the additional information allows one to test particular hypotheses regarding haemorheology (see Sect. 3).Interestingly, the seminal paper by Lipowsky et al. [15] already hinted at this: they observed fairly large discrepancies when comparing the predicted pressure drops with the values obtained experimentally.In some regions, pressure drops were up to six times greater than expected.They suggested that the non-Newtonian behaviour of blood was the main reason for the discrepancy.With the combined availability of both topology and flow, these issues can be investigated at last.

Graph representation
The most convenient method to describe networks is in the form of a graph.Several previous studies already explain clearly how this concept can be applied [6,11,23], so only a brief description is given here; the focus is on the differences in the solution methods if also flow information is available through experiments.

A simple example network
To introduce the nomenclature and tools for graphs, we here describe a very simple (vascular) network, as shown schematically in Fig. 3.The network contains a number of nodes (also referred to as vertices), connected by branches (also referred to as edges).The nodes and branches are identified by an index, i ¼ 1. ..n and j ¼ 1. ..m, respectively; in this simple example n ¼ 5 and m ¼ 6.The pressure is defined at each node,1 denoted as an element of vector P i .
To described the topology of the network, i.e. which nodes are connected by which branches, several options exist.The most straightforward is the use of a connectivity matrix, C ij .This (n Â m) matrix contains elements c ij , with i the node and j the branch.If a particular node i is not connected to a particular branch j, the value of the element c ij is zero.When the direction of the flow (i.e.into or out of nodes) is known, the element values are either -1 or 1.By convention, a negative value denotes that flow leaves the node.Such a connectivity matrix then describes a so-called directed graph for a given flow network [11].If the flow direction is not yet known, the undirected connectivity matrix only contains ones and zeroes; the direction of the flow will then appear later on in the analysis (e.g in the sign of the flow rate through a particular branch).For the simple model network, the directed connectivity (C d ij ) matrix is given as: Note that the connectivity matrix is generally sparse, as nodes are usually only connected to three of the possibly many branches.This sparsity is an important characteristic for the practical solution implementation when large networks are considered.
The flow rate in each branch will be denoted by the (1 Â m) column vector Q j .For our simple example network, it has the following elements: As blood is an incompressible fluid, there can be no accumulation of mass (or volume) at each node.The sum of all flow rates entering and leaving a node must Fig. 3 A simple network, used to illustrate the nomenclature of graphs thus be zero.To facilitate bookkeeping, we can define the source vector S i , with dimension (n Â 1) as S i ¼ C d ij Q j .For our example: For nodes i that are not at the edge of a network S i should be zero to ensure conservation of mass, e.g. for node 2 in Fig. 3: The flow Q j will be the result of a pressure difference between the two connected nodes, DP j .Again, we can use the connectivity matrix to find an expression for this pressure difference from the pressure at the nodes: As example, DP 2 ¼ P 3 À P 2 in the simple example network.Note that the hypothesis that flow is a consequence of a difference in pressure, i.e. a gradient of a potential, ensures that Kirchhoff's second rule (in our analogy: no pressure drop in a closed loop) is satisfied.
To link the flow rates in the branches with the pressure in the nodes, a model has to be constructed that describes the resistance for each branch.This resistance will likely be a function of the length and diameter of the branch, the flow rate and the rheology.In Sect.2.5 a more detailed discussion is given; here it is assumed that the resistance of a branch (R j ) can be calculated based on Poiseuille's law: For the pressure drop due to a given flow, we thus have: Note that here Hadamard (element-wise) multiplication is implied, so all three terms are vectors of the same size.From Eqs. 4 and 6, we have: To find the pressure at each node for a given network with known flow rates (and thus also resistances), we can solve this system of linear equations for P i .As C d ij is generally not square, a solution cannot be found using its inverse.For the general case, the Moore-Penrose pseudo-inverse can be used.The presence of measurement errors requires the use of an optimisation process (e.g. using least squares) to find an estimate Pi for the generally overdetermined system of equations.Numerical solution of the system of equations, as well as all data processing, is performed in MATLAB (R2013b, The Mathworks) in this study.
If the flow rates are not known, as is the case for most of the previous studies reported in Sect.2, a solution can still be found due to the fact that conservation of mass needs to hold, S i ¼ C d ij Q j ¼ 0. This equation can be combined with Eq. 7 to find both P and Q, see e.g.Kassab et al. [11].

Boundary conditions
Equation 7 represents the system of linear equations that describe the pressure differences between the nodes.Unless at least one reference pressure is known at a given node, all pressure will be relative to an unknown reference pressure.If known, this reference pressure can easily be prescribed by appending the value to the right-hand side vector and adding a row to the connectivity matrix with a single non-zero element at the appropriate node.
When pressure and flow need to be found simultaneously, missing or invalid branches will significantly alter the end result.However, here the flow rates are available from the experimental data.Therefore such invalid (or missing) segments are simply ignored (or absent) in the solution procedure.For instance, in Fig. 3 any segment, apart from j ¼ 1, can be removed and the network will remain connected and the pressure distribution can be determined correctly.This also means that 'dead ends' -at the edge of the network or elsewhere-do not require special treatment.Naturally, when more information is available in the overdetermined system, the effects of measurement uncertainties will be reduced.

Building the network
For simple data sets containing only a few connected blood vessels, the graph is best constructed manually.
Naturally, this is no longer feasible in larger data sets.Fully automated image data segmentation to obtain a graph representation has been described in myriad previous studies [9,23,26].Note that again these methods are generally based only image data (c.f. the left hand image in Fig. 2).However, the velocity data (such as the right hand figure in Fig. 2) can facilitate the detection of vessel segments.A detailed example of this approach is given in Kloosterman et al. [13]; here only the key points are summarised in the following paragraph.
In Fig. 4, an example is given of segmentation based on velocity field data.In the top left figure, the input is shown: a time-averaged velocity vector field typically containing hundreds blood vessel segments.The figure on the top right shows a detailed view of the vector field, with the velocity magnitude shown as grey-scale background.As a first step in the processing, a mask is created by thresholding the velocity magnitude, using a threshold value comparable to the measurement error.This mask is a binary matrix with the same size and resolution as the velocity field; all image processing steps use this resolution.This mask is shown as the outline in the bottom figures.The socalled skeleton of this mask is constructed, shown as the colour-coded line segments in the bottom left figure.At locations where the skeleton/centerline bifurcates a node is defined, here denoted by an open circle.Nodes and the segments in between nodes are then numbered (the colour coding in the bottom right represents the segment index).The result of these steps is a collection of nodes and segments.An extensive set of these collections for various embryos at different developmental stages is available in the aforementioned paper [13].This data will serve as the starting point for the further processing in this study.
The relationship between the (arbitrarily-numbered) segments and nodes are captured in the connectivity matrix.For each node, the segments are found that connect to this branching point (note the three non-zero entries around each node in the colourcoded skeleton image shown in the background of the bottom right of Fig. 4).As the velocity field is also available, we can directly constructed the directed connectivity matrix, C d ij .To find the flow direction, the velocity along a few points on the skeleton is evaluated in the original vector field: if this flow is on average toward the node, the value c ij is 1; flow away from the node results in a value of À1.This process is repeated for all nodes to construct the connectivity matrix.A graphical representation of a small part is shown in Fig. 5.Note that here the c ij values have been multiplied with the flow rate through the branch, Q j (see later).Red denotes positive flow, i.e. toward a node, while blue denotes flow away from a node.The sum of each row in this matrix should be zero, see e.g. the elements indicated by the rectangle that represent the situation of node 35: branch 39 enters the node, while branch 63 and 64 leave it.End-points can also be defined, these are nodes that are only connected to a single branch.Identification of these end-nodes is important if particular boundary conditions are specified.
While this automated processing is fairly robust, occasionally 'gaps' in the network can occur.This can be due to data drop out (vessels may not appear continuous in the source data) or by erroneous assignment of nodes (e.g. two closely placed nodes, instead of the correct single node at a bifurcation).This may lead to two unconnected graphs, which corresponds to an underdetermined system for Eq. 7.This can easily be fixed manually after visual inspection of the graph representation.

The hydraulic resistance
With the connectivity matrix described above, one half of the graph is available.The other half entails assigning a hydraulic resistance to each of the vessel segments or branches.As stated earlier, this resistance is a lumped parameter describing the total pressure drop in the fluid going from one node to another.The exact value will be dependent on the geometry (of the branch, but also the junctions at its ends), the conditions of the flow and the rheology of the fluid.
In Eqs. 5 and 6, the pressure drop is estimated based on Poiseuille's law.This may seem like an overly simplified approach.However, it turns out that the flow condition in the vasculature under investigation satisfies the criteria for Poiseuille flow.The velocity rarely exceeds 0.5 mm/s and vessels are generally smaller than 0.2 mm [13].With an approximate kinematic viscosity of 2-3 mm 2 =s, the Reynolds number (Re UD=m) is much smaller than unity in all branches.This indicates that viscous forces will dominate over inertial terms.This implies that many inertia-related phenomena (e.g.due to pulsatile flow, vessel curvature or entrance effects) are negligible.For a detailed discussion, including the appropriate dimensionless numbers, one is referred to [19].
As an illustration, we here evaluate the additional energy dissipation due to entrance effects, which are expected to be the most prominent additional losses in the present case.The losses occur due to the fact that a flow requires some distance after a geometry change to approach the fully-developed parabolic profile.Expressed as fraction of the Poiseuille pressure drop, the additional loss can be estimated using the dissipation ratio Z entry [17]: The diameter over length ratio of the vessel, D / L, is generally smaller than unity.The average Reynolds number in the network as shown in Fig. 4 is of order 0.01.Therefore, the additional pressure drop is thus only a small fraction of the Poiseuille pressure drop, so we can safely ignore entrance effects here.For other applications, in particular in larger vessels with higher Reynolds numbers, they can be incorporated using a correction based on Eq. 8.
Implicitly, it has been assumed that the rheology of blood can be captured by a single parameter: the viscosity (l).While human blood is known to exhibit shear-thinning behaviour, it should be stressed that the current data set is obtained in a chicken model system.Avian blood, in particular in the embryonic stage, is very close to a Newtonian fluid [1,19].In Sect.3.2, the possibilities of having a different value of the viscosity for each branch is explored.
Based on these considerations, it can be concluded that Poiseuille's law is here indeed appropriate to capture the pressure drop in a vessel branch.Furthermore, the average velocity field data contains sufficient information to model each branch.The relevant parameters for Poiseuille's law-diameter, flow rate, length-are extracted branch by branch: the velocity profile is described using a parabolic fit to the data, perpendicular to the branch skeleton.This is repeated along the downstream direction and the results are averaged to find a robust estimate of the diameter D j and mean velocity V j;mean for each branch. 2he flow rate is then calculated from these two parameters (Q j ¼ p 4 D 2 j Á V j;mean ).Lastly, the vessel length along the skeleton (or centerline), L j is stored, which is approximately 10% longer than the Euclidian distance between the end points [13].This procedure provides all terms required in Eqs. 5 and 6 for each branch, so that the pressure distribution in the network can be found from Eq. 7; note that the latter system of equations is linear due to the use of Poiseuille's law, as pointed out by Pries et al. [22].

Two application examples
To illustrate the capabilities of the tools introduced in the first part of this manuscript, two examples are given here.They both focus on utilizing the additional information that is available, i.e. flow data instead of just imaging/topology.

Iterative correction of divergence
Experiments such as shown in Figs. 2 or 4 inevitably contain measurement noise.More specifically, the velocity data will contain an error.This error will propagate through the data processing steps, so that the graph network as shown in the bottom half of Fig. 4 will contain errors in the values for Q j and D j (which are both derived from a fit to the velocity profile).Errors in the connectivity, as described in the end of Sect.2.4, are assumed to be absent.
The network approach can assist in identifying and correcting the errors in the branch properties.This is particularly useful because for this type of in vivo measurements it is unfeasible to perform reference or control measurements [18].The main idea behind the identification and correction process is the evaluation of the mass balance at each node.As defined in Sect.2, S i should be zero for all nodes that are not end nodes.
The error at individual nodes can be visualised, as shown in the left hand side of Fig. 6.In this figure, the discs show the absolute value of the sink term at each node, i.e. a larger disc denotes that more fluid is unaccounted for.The total amount of fluid unaccounted for in the entire network shown in Fig. 6, RjS i j, is 14:4 Â 10 À3 mm 3 =s, with the largest sink being 7:6 Â 10 À4 mm 3 =s (the node in the centre, near x = 1.6 mm, y= 2.8 mm).To put these errors in perspective, the total flow through the network, RQ j , is 0:14 mm 3 =s-only an order of magnitude larger.The flow rate averaged over all branches is 4:9 Â 10 À4 mm 3 =s.Visual inspection of the data near the largest sink showed that the cause was a branch with a considerable velocity, yet a relatively small diameter; the automated processing underestimated the flow rate in this particular branch.The analysis can be done more refined, e.g. by looking at the relative local flow balance, ðQ out À Q in Þ=Q in .Furthermore, these errors can also be evaluated for different branch diameter groups, to further study the underlying causes of the errors.The ability to exclude branches that are very likely erroneous will improve statistical analysis in physiological studies.
In the previous paragraph nodes were evaluated individually to assess the accuracy of the flow in the branch.However, in the present data set and framework it is possible to evaluate nodes and branches in their context, rather than as isolated entities.To highlight this, an iterative scheme is here proposed to reduce the measurement errors in the flow rates.To illustrate the scheme, consider branch 2 in the simple network shown in Fig. 3.This branch is connected to two nodes (2 and 3).At each node conservation of mass holds, so we can write: We can find an alternative value for Q 2 by finding the flow rate that would perfectly satisfy conservation of mass at the two nodes: We can choose to either completely replace the old value with the new estimate or use some sort of overrelaxation approach to update its value incrementally (e.g.Q j !kQ j;new þ ð1 À kÞQ j ).This procedure is referred to as Successive Overrelaxation [24].An alternative visual explanation of the approach can also be found using a real connectivity matrix, as shown in Fig. 5.The values along the rows within the rectangle should add up to zero.We can re-evaluate the second value (branch 63), but this will also mean that another row is affected, as this same branch also occurs there, with an opposite sign (as indicated by the vertical dashed line).
To process an entire network, the following algorithm can be used: 1. Randomise order of list of branch indices, j ¼ 1. ..m

2.
Select a branch j

3.
Find the two nodes that are connected to this branch

4.
Find the flow rate that satisfies conservation of mass at both nodes

5.
Update the flow rate, Q j !Q j;new

6.
Repeat from step 2 until all branches have been updated 7. Repeat from step 1 until the sum of sinks RjS i j is converged The main motivation behind the scheme is that the random error in the four branches cancel out to some extent compared to the error in the original flow rate.In a way, it can be interpreted as a smoothing process, but it will only smooth information that is nonphysical, i.e. divergence in the flow field.Care must be taken when there are exceptionally large values (e.g.flow rates that are orders of magnitude too high).Such overestimates will be redistributed over the network, increasing the overall flow rates.
To demonstrate the effectiveness, the network shown in the left-hand side of Fig. 6 was processed using the proposed iterative correction scheme.The result is shown on the right-hand side of the same figure.The discs at each node again show the magnitude of the sink term.As can be seen in the figures, the sinks are greatly reduced.The colourcoding of the branches indicate the change in each branch: blue denotes an increase in flow, red denotes a decrease, and white branches are not changed at all.The thickness of each branch indicates the flow rate in both cases.Various tests using different starting branches and relaxation values (k) appeared to converge to the same solution.The results shown here were obtained for k ¼ 0:5 and 200 iterations.Computational efforts were minimal (less than a minute on a desktop PC) even for a naive implementation of the algorithm.For more complex networks and models, i.e. non-Poiseuille resistance terms or coupled flow / mass-transfer models it will be useful to optimise the numerical solution procedure [24].It cannot be expected that iteration will always converge to the 'true' solution, so visual inspection is recommended.Changes to this particular network were relatively minor, with only a handful of branches required significant flow correction.
To quantify the improvement, the total sum of the magnitude of the sink term is shown in Fig. 7, together with the results for an alternative data set.The dashed line represents the data set shown in Fig. 6, while the other data set was taken in the same embryo and region, but at a different stage during development.The total sink term, RjS i j, reduces from 0.0144 to 0.0021 mm 3 =s, a reduction of 85 %.For the other data set a similar decrease is observed.Note that comparison of the absolute values of the two cases is difficult, because not only the imaging conditions-and thus measurement error-are different, but also the number and type of blood vessels in the field of view.
The iterative scheme proposed shows a great improvement in the network in terms of satisfying the conservation of mass principle, so it is expected that the corrected networks better represent the real situation.Naturally, this proposed scheme is the most simple approach possible and several refinements can be introduced.For instance, additional restrictions can be imposed, e.g. if a reference value for a particular branch is known or to ensure that the total flow through the network may not change.Alternatively, one can only correct branches that are 'suspicious', i.e. branches that have two nodes with significantly strong sink terms.
As a final note on this scheme, it must be mentioned that the values of Q j have been corrected.These values were determined from the velocity field, via a mean velocity and diameter in each branch.The updated values for Q j can in turn be used to correct V j;mean or D j if needed, for instance if wall shear stresses are studied.Which of the two should be corrected depends on the specifics of the method by which the data is obtained and processed.

Haemorheological model testing
The previous example was based on Kirchhoff's first law, here implying conservation of mass at nodes.The second application example to demonstrate the unique capabilities of network flow data is based on Kirchhoff's second law.Using Eq. 7 the pressure distribution in a network can be found (see also e.g.Fig. 9, discussed in more detail later on).This has been reported before in earlier studies, using topology data only.However, as in this case the flow rates are already available, they do not have to be approximated in the solution procedure too.This means that one layer of uncertainty is removed in the resulting pressure field.Furthermore, the availability of flow rates provides additional opportunities, as discussed below.
In Sect.2.5, it was assumed that the viscosity was constant throughout the network.However, it has long been established that the effective viscosity in a branch can be highly variable.The underlying causes are mostly related to the spatial distribution of erythrocytes.Their volume fraction (and thus the viscosity) varies per branch-the Fa ˚hraeus effect [22,23].Furthermore, the presence of a cell-free layer reduces the effective viscosity, in particular in smaller branches, the Fa ˚hraeus-Lindqvist effect.Various models have been proposed that incorporate these effects, but there is a scarcity of data to properly validate them.Furthermore, it is expected that there is a wide variation between species and even during development [19], so it is not clear if haemorheological models can be used across species.In the following, a method is described that may not be able to be used to determine rheological behaviour from scratch, but it can help in comparing various models.
To illustrate the approach, we once more consider the simple network shown in Fig. 3. Node 2 and 4 are connected by two routes: via branch 4 and via branches 2 and 3.The pressure difference between the two nodes, P 2 À P 4 , must be path-independent.Using Poiseuille's law, the pressure difference can be expressed as follows: The resistances of branch 2 and 3 could be added as they are resistances placed in series.As stated, the pressure difference between both routes should be identical.The difference can therefore be used as penalty function: As L j , D j and Q j are all known, this leaves only the viscosity in each branch unknown.Even in the hypothetical case where DP 0 is exactly zero (no measurement noise), the absolute value of the viscosity cannot be established from Eq. 12.To do this, a known pressure drop is required as additional information.However, it is possible to test various haemorheological models this way.By substituting a particular model, l j ¼ f ðD j Þ, the penalty function (Eq.12) can be evaluated.The function with the minimum penalty function will best describe the Fig. 7 The total sum of the sink terms, RjS i j, during the iteration process for two different data sets rheology.As stated, the absolute value cannot be determined, but only the relative dependency on D.
To utilise the entire network (instead of the single penalty function shown in Eq. 12), the following approach is proposed: The pressure at each node is first estimated from Eq. 7. The values R j are obtained using the known values of L j , Q j , D j and the viscosity l j using the model under investigation.Pi is estimated using the pseudo-inverse method and by specifying a reference pressure at one node.
The difference between the prediction based on the 'Poiseuille' estimate and the approximation P represents the performance of the rheological model: The penalty term DP 0 has here been redefined to denote the difference between the prediction 'a priori' based on Poiseuille's law for a particular branch and the result of the matrix inversion approximation for that same branch.Note that the estimated values for the pressure ( P) are now used in the first term of the right-hand side.To quantify the total discrepancy, the standard deviation of all values of DP 0 is calculated.This value is normalised using the mean pressure drop.This normalised error will be denoted E.
To demonstrate the use of this approach, two rheological models are compared: (1) a constant viscosity and (2) a model incorporating the Fa ˚hraeus-Lindqvist effect.For the latter we make use of the model given by Pries et al. [22], who provide a fit based on a large number of experiments: In Fig. 9 the resulting pressure distribution is shown for the vascular network.The reference pressure was set to zero at the node in the main branch that enters the field of view (P 0 ¼ 0).For this particular case the diameter-dependent rheology model by Pries et al. was used.The total pressure drop over the fieldof-view is approximately 3 Pa; for reference, the mean gauge pressure in a vitelline artery, i.e. the largest branch in this type of vasculature, at this stage (HH14, around 50 h of development) is 50 Pa [8].Only a small fraction of the total available driving pressure is thus dissipated in the modelled network.
The normalised error E is 0.461 for the constant viscosity case, while it is 0.408 for the more complex model -a decrease of 12 %.Note that this decrease is not due to the slightly lower mean viscosity, as this is taken into account by the normalisation with the mean pressure drop.The lower value for the diameterdependent model suggests that it better describes the true situation.This very simple comparison is by no means presented as a serious haemorheological study, it just serves to illustrate the procedure.

Concluding remarks
In this manuscript analysis tools are presented for the data that has become available with the latest Fig. 8 The two haemorheological models: (1) constant viscosity and (2) diameter-dependent viscosity as described by Pries et al. [22].Also shown, in the bottom part, is a histogram of the branch diameters, based on a total of 271 vessel segments generation of measurement modalities, in particular the ones that present both topology and flow information.The tool set, based on graph representation, opens up unique opportunities, as blood vessels can not only be studied as separate entities, but in their context.An important feature is the significant reduction in data, from large velocity fields to a compact representation using matrices.This reduction, together with straightforward matrix manipulations, facilitates the study of the haemodynamics in large networks.Two concrete examples are presented that make use of these unique features: (1) the identification and correction of measurement errors to ensure conservation of mass and (2) the ability to test various haemorheological models via the pressure field.For both examples a proof-of-princple result is shown using data obtained in a chicken embryo model system; it can readily be applied to other data sets.

Limitations and outlook
Apart from presenting the opportunities, it is also relevant to evaluate the limitations of the network approach.These may arise from the imperfections in the input data, but also from the assumptions that were made in the model formulation.For the former, it is obvious that the experimental data will not always perfectly reflect reality.There will be random measurement noise in the velocity fields and thus also in the derived branch parameters.In particular the uncertainty in the vessel diameter will have a major impact on e.g. the pressure distribution, as the pressure drop is proportional to the diameter to the fourth power.Furthermore, there can be systematic errors due to three-dimensionality of the network or differences in orientation between the measurement plane and the center-plane of the network.The will give rise to unphysical results: for instance, the flow rate in a branch may change along the downstream distance or the flow balance at each node may no longer be zero.If this is the case, the measurements and model need to be extended to three dimensions.
For the microrheological testing procedure described in this manuscript, a major limitation is the fact that it is a 'trial-and-error' approach.It cannot be used to determine rheological behaviour, but only to evaluate the (relative) performance of existing models.This is due to the lack of reference pressure drop measurements.The small scale, intrusive nature and relative low pressure drops make such measurements far from trivial.
The simplicity of the network description is the result of a series of assumptions about the flow in the vasculature.If the present case, it has been argued that Poiseuille's law is valid.In many other cases this will no longer be the case.For instance, in capillaries the behaviour of individual formed elements must be taken into account and the fluid is no longer a homogeneous medium.As another example, the assumption that entrance effects are negligible will not be valid for networks with higher Reynolds numbers.The same holds for networks observed during very early stages, where it can be difficult to identify distinct individual branches [14].A porous media approach might be a better alternative here.
Fortunately, most of these limitations can readily be addressed by an appropriate extension of the model, e.g. by incorporating entrance effects using Eq. 8, by Fig. 9 Pressure distribution in the network, using the diameterdependent haemorheological model.The pressures are relative to the reference pressure P 0 ¼ 0 specified at the bottom of the network incorporating local hematocrit variation in the capillaries [22] or formulating the approach for threedimensional structures [26].
The network method presented here can further be extended, for instance to look at mass transfer or transit times in the network.Future work will first focus on the systematic analysis of the large experimental data sets that are available [13], in particular to see whether rules can be formulated that decide the fate of particular vessels.This will hopefully shed light on the intricate processes that guide the flowinduced modifications in a network during cardiovascular development.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/),which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Fig. 2
Fig. 2 An example of a state-of-the-art measurement of the flow in vascular networks.(Left) Raw brightfield image data series.(Middle) Mean velocity field obtained by particle image

Fig. 4
Fig. 4 Reduction of the vector field to a graph representation; the right hand figure shows a close-up of a section of the total field-of-view.The vectors in the bottom right denote the general

Fig. 5 A
Fig. 5 A graphical representation of the connectivity matrix of the vessel network shown in Fig. 4. Some elements of the matrix C d ij Á Q j are shown, with the colour-coding signifying flow into or out of the node (red and blue, respectively).(Color figure online)

Fig. 6
Fig. 6 Network flow consistency before and after iterative correction.The discs mark the node locations, with the size representing the magnitude of the sink terms.The line segments represent the branches; colour-coding represents the change in with D in lm.Both models are shown in Fig.8.In this figure, the dashed line represent the constant viscosity model (here l = 3.2 mPa s is chosen).The model by Pries et al. approaches the same value for large D. For smaller values of D, there is a minimum in the viscosity around 7 lm; for smaller values the viscosity increases dramatically, as here the blood vessel diameter becomes smaller than the typical size of an erythrocyte.Also shown is a histogram of the branch diameters (based on a total of 271 branches and a bin size of 20 lm).The dots indicate the value of the viscosity for each branch.As can be seen, for the diameters under consideration in this network the viscosity varies by nearly 20 %, with all branches larger than the minimum-viscosity diameter.