Computer simulations reveal complex
distribution of haemodynamic forces in a
mouse retina model of angiogenesis
By Miguel O. Bernabeu1,2 , Martin Jones3 , Jens H. Nielsen4 , Timm
Krüger5,2 , Rupert W. Nash2 , Derek Groen2 , Sebastian Schmieschek2 ,
James Hetherington4 , Holger Gerhardt3 , Claudio A. Franco6,3,∗ ,
Peter V. Coveney2,∗

arXiv:1311.1640v2 [] 15 Jul 2014

1

CoMPLEX, University College London, Physics Building, Gower Street, London,
WC1E 6BT, UK
2
Centre for Computational Science, Department of Chemistry, University College
London, 20 Gordon Street, London, WC1H 0AJ, UK
3
Vascular Biology Laboratory, London Research Institute - Cancer Research UK,
Lincoln’s Inn Laboratories, 44 Lincoln’s Inn Fields, London WC2A 3LY, UK
4
Research Software Development Team, Research Computing and Facilitating
Services, University College London, Podium Building - 1st Floor, Gower Street,
London, WC1E 6BT, UK
5
Institute for Materials and Processes, University of Edinburgh, King’s Buildings,
Mayfield Road, Edinburgh EH9 3JL, Scotland, UK
6
Instituto de Medicina Molecular, Faculdade de Medicina Universidade de Lisboa,
Lisboa 1649-028, Portugal
∗
Equally contributing senior authors

There is currently limited understanding of the role played by haemodynamic forces
on the processes governing vascular development. One of many obstacles to be overcome is being able to measure those forces, at the required resolution level, on vessels
only a few micrometres thick. In the current paper, we present an in silico method
for the computation of the haemodynamic forces experienced by murine retinal
vasculature (a widely used vascular development animal model) beyond what is
measurable experimentally. Our results show that it is possible to reconstruct highresolution three-dimensional geometrical models directly from samples of retinal
vasculature and that the lattice-Boltzmann algorithm can be used to obtain accurate estimates of the haemodynamics in these domains. We generate flow models
from samples obtained at postnatal days (P) 5 and 6. Our simulations show important differences between the flow patterns recovered in both cases, including
observations of regression occurring in areas where wall shear stress gradients exist.
We propose two possible mechanisms to account for the observed increase in velocity and wall shear stress between P5 and P6: i) the measured reduction in typical
vessel diameter between both time points, ii) the reduction in network density triggered by the pruning process. The methodology developed herein is applicable to
other biomedical domains where microvasculature can be imaged but experimental
flow measurements are unavailable or difficult to obtain.
Keywords: angiogenesis, mouse, retina, blood flow, simulation, shear stress,
image-based, lattice-Boltzmann
Article submitted to Royal Society

TEX Paper

2

Bernabeu, Jones, Nielsen, et al.

Figure 1. Murine retinal vascular plexus 6 and 21 days postnatal (left and right panel,
respectively). Within days, the primitive vessel network remodels into mature vasculature.
Samples were collected, mounted, and imaged as described in Section 3.

1. Introduction
Despite recent advances in vascular biology, the mechanisms underpinning vascular
development remain poorly understood. It is therefore crucial to gain further insight
into the mechanisms governing the formation of complex vascular networks and
their response to external stimuli. The translation of these results holds the key to
the improvement of therapies modulating vascular patterning and sprouting for the
treatment of stroke, ischaemia, retinopathies or cancer, the leading cause of death
worldwide.
One of the pressing questions in the field is establishing how primitive vessel
networks remodel into a hierarchically branched and functionally perfused network
of arteries, arterioles, capillaries, and venules (Figure 1). In recent years, the main
molecular mechanisms regulating endothelial cell behaviour during vessel formation
have been elucidated using experimental techniques (Jones et al., 2006; Potente
et al., 2011). However, important challenges remain: i) understanding how cell-level
mechanisms integrate to give rise to systems-level behaviour and, ii) understanding the impact in vascular patterning of the interplay between cellular molecular
regulation and haemodynamic forces (i.e. vascular mechanotransduction). These
problems are hard to address due to the multiscale and multiphysics nature of the
processes involved. Systems-level behaviour arises from highly non-linear, tightly
coupled interactions between subprocesses at different spatial and temporal scales.
Furthermore, it has been recently proposed (Bentley et al., 2013) that a tighter
integration between experimental and computational work is required in order to
tackle these questions. Working in a feedback loop, computational models should
be capable of generating new hypotheses, rather than merely reproducing experimental data. In turn, experiments should provide new biological insights based on
these hypotheses and help to further refine computational models.
Multiple animal models have been proposed for the study of vascular development. Examples include the mouse retinal and embryonic vasculature (Geudens and
Gerhardt, 2011), zebrafish vasculature (Lawson and Weinstein, 2002; Siekmann and
Lawson, 2007), and hyaloid vasculature (Lobov et al., 2005). In recent years, there
has been increasing interest in the development of in silico models for the close inArticle submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

3

spection of certain vascular developmental aspects. To date, most work concerning
simulation of retinal haemodynamics for the study of vascular mechanotransduction
(see Section 2c for a review) has suffered from a number of limitations including: i)
limited availability of spatial information due to the use of low resolution imaging
modalities, ii) oversimplification of the haemodynamics by considering the retinal
plexus to be a network of one-dimensional vessel segments, and iii) unavailability
of the computer code developed. We believe that the model simplifications cited,
although appropriate in some applications, may fail to capture complex flow patterns important for understanding the interplay between molecular regulation and
haemodynamics during development. Hence, in the current work we introduce a
computational workflow — and make the source code available — aimed at generating in silico estimates of the haemodynamic forces acting on samples of mouse
retinal vasculature imaged during development, typically within the first postnatal week. The workflow involves the following steps. First, high-resolution scanning
confocal microscope images are obtained and segmented in order to generate a
binary mask of the vessel lumen. Second, luminal centrelines and radii are computed in a process known as skeletonisation. Next, three-dimensional (3D) models
of the luminal surface are reconstructed based on the computed skeleton. Finally
blood flow simulations are run in order to obtain estimates of blood velocity and
wall shear stress with an open source highly-parallel computational fluid dynamics
(CFD) solver, known as HemeLB (Mazzeo and Coveney, 2008).
The purpose of this paper is therefore three-fold. First, to describe the computational methods developed and to survey the literature for data not accessible in
our experiments but necessary for model set-up. Second, to validate our methods in
simplified scenarios where analytical solutions are known. Third, to present and to
analyse simulations in order to gain insight into the dynamics of retinal blood flow
during development. The paper is structured as follows. In Section 2, we survey the
literature for previously proposed models of retinal flow and for the experimental
data necessary to set up our simulations. Next, in Section 3 we present the methods
used for image processing and 3D model reconstruction as well as the validation
methodology adopted. Section 4 presents the main results on model reconstruction,
validation, and a set of simulations on the reconstructed 3D models. Finally, Section 5 summarises the main contributions of the work and outlines the areas where
we plan to apply the computational pipeline developed.

2. Retinal vascular structure and flow
(a) Vascular structure and its development
Angiogenesis defines the formation of new blood vessels from pre-existing ones
and can be split into two distinct phases: sprouting and remodelling. During sprouting, new vessels form and invade avascular ischaemic areas, where tissues experience
hypoxia and nutrient deprivation. This process is modulated by the secretion of various growth factors, including vascular endothelial growth factor (VEGF), through a
cascade of signalling events. The end point of this phase is the formation of a highly
branched and poorly perfused network of capillary connections. Remodelling is responsible for the creation of a hierarchically branched and efficient vascular tree,
containing defined arteries and veins and an optimised vascular capillary network. A
Article submitted to Royal Society

4

Bernabeu, Jones, Nielsen, et al.

vital step during vascular remodelling is the removal of redundant vessel segments:
vessel pruning. Importantly, angiogenesis is a very dynamic process occurring not
only during development but also in adulthood (e.g. wound healing and tumour
formation).
The neonatal mouse retina has become one of the main experimental models
for the study of the mechanisms involved in blood vessel development and patterning (Fruttiger, 2007; Uemura et al., 2006; Gariano and Gardner, 2005). The mouse
retina is avascular at birth and develops through a consistent series of events. Astrocytice (a type of glial cell) and neuron derived vascular endothelial growth factor A
(VEGFA) stimulates sprouting angiogenesis from pre-existing blood vessels at the
optic nerve. Under a gradient of VEGFA, vessels expand radially in the superficial
layer of the retina, in a very characteristic pattern (Potente et al., 2011). It takes
this process around 8 days to cover the entire surface of the mouse retina. Vascularisation of the superficial layer is followed by a second phase of sprouting, where
endothelial cells from the superficial venous plexus sprout and penetrate the deeper
layers of retina to form, firstly, a deep and, secondly, an intermediate capillary bed
(Brown et al., 2005). The vascular plexus finally matures about 20 days after birth.
In the mature mouse retina, vessels are found to be predominantly arteriolar in
the superficial layer and predominantly venular in the deep capillary bed (Paques,
2003). The artery feeding the retina arrives at the optic disc and divides into 8 to
9 radiating retinal arteries. These arteries, with luminal diameter of up to 28 µm
(Ninomiya and Inomata, 2006), side branch into smaller arterioles at close to 90
degrees angles from the parent arteries. The arterioles (10 µm to 12 µm in diameter
(Ninomiya and Inomata, 2006)) take a relative long course before abruptly changing
direction to run towards the intermediate and subsequently deep capillary beds
(5 µm to 6 µm in diameter (Ninomiya and Inomata, 2006)). Some authors (e.g.
Paques (2003)) have suggested that the superficial layer is mostly capillary-free
with only a few direct connections between arteries and veins.

(b) Haemodynamics
The relationship between haemodynamics and pathogenesis of various eye disorders has prompted researchers to analyse retinal blood flow in both basic and
clinical research domains. Early examples are the work of Feke et al. (1989), who
measured total retinal blood flow and its regional distribution in humans, and Alm
and Bill (1973) who studied blood flow rates in various tissues of the primate eye.
Later advances in imaging techniques (e.g. optical coherence tomography and related modalities, laser Doppler velocimetry) have expanded our understanding of
retinal haemodynamics. High resolution in vivo measurements of retinal flow have
been obtained in various species: mouse (Brown et al., 2005), rat (Zhi et al., 2011),
and human (White et al., 2003; Wang et al., 2009). Several authors (White et al.,
2003; Yazdanfar, 2003) presented evidence of the pulsatile nature of retinal blood
flow despite early claims (Alm and Bill, 1973) that only retinal arteries — and
not veins — exhibit systolic to diastolic flow rate variations. Further developments
enabled quantitative analysis of typical arterial and venous flow in both healthy
(Wang et al., 2007) and diseased (Wang et al., 2009) human retinas as well as
during development in mice (Brown et al., 2005). Table 1 compiles some of the
Article submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

5

measurements of murine and human retinal vessel diameter, velocity, and flow rate
available in the literature.
Of relevance to our study is the work by Brown et al. (2005), who obtained in
vivo measurements of blood flow in various parts of the murine eye including the
retina from birth to postnatal day (P) 16. Table 2 presents some of their findings.
A clear trend of increase in retinal blood flow after P3 is observed. The authors
attribute the large variability of the results (note standard deviation in Table 2) to
the natural variation in the time course of the remodelling processes involved.
Several authors have also studied the typical pressure difference driving flow
in the retina, the so-called ocular perfusion pressure (OPP). The pressure at the
central retinal artery is often approximated with mean arterial pressure (MAP)
measurements at eye level (e.g. carotid arterial pressure (Wright and Harris, 2008),
subclavian artery (Hardy et al., 1994)). The effective venous pressure is considered
equivalent to the intraocular pressure (IOP) (Kiel, 2010). Table 3 summarises the
values of MAP and IOP reported by several authors. OPP values of ∼57 mmHg
are consistently reported across species. Finally, retinal blood flow is known to be
autoregulated by the modulation of retinal vessel compliance in response to MAP
and IOP variations (Paques, 2003; Hardy et al., 1994).
From a rheological point of view, blood is a shear-thinning fluid (i.e. its viscosity
is a decreasing function of shear rate (Chien, 1970)). When flowing at sufficiently
large shear rates (typically greater than ∼ 1000 s−1 ) blood can be modelled as a
Newtonian fluid (i.e. constant viscosity) with no significant effect on the simulated
haemodynamics (see Bernabeu et al. (2013) and references therein). However, at
lower shear rates, viscosity quickly increases due to e.g. red blood cell (RBC) aggregation. Our work concerns simulation of blood flood in small arterioles, venules,
and capillaries where shear rate is expected to be lower than the aforementioned
threshold. In fact, Nagaoka and Yoshida (2006) measured shear rates as low as
(606 ± 115) s−1 in the venules of the human retina. Therefore, we will take the
shear-thinning properties of blood into account in order to improve the fidelity of
the shear stress computed in our model. Table 4 compiles some of the values of
murine blood viscosity as a function of shear rate available in the literature. Other
rheological properties derived from the presence of RBCs, such as the FåhræusLindqvist effect (see e.g. (Pries et al., 1996)), will not be considered.
Finally, we note that experimental data on haemorheology changes during development is limited. Windberger et al. (2005) observed a steady increase in blood
viscosity in rabbits and cats from fetal stages to adulthood. The changes were
more pronounced within lower shear rates regimes: from 3.00 mPa s to 9.29 mPa s
at 0.7 s−1 and from 2.48 mPa s to 3.62 mPa s at 94 s−1 during the first 30 days of
life in rabbits. Due to the scarcity of available data, we will not include this effect
in our model.
The experimental data on pressure distributions and haemorheological properties surveyed in this section will be used to set up our flow simulations. This is
done due to the impossibility of obtaining such information directly from our experimental model. The data summarised in Table 1 will be used to validate our
experimental measurements of vessel diameter and in silico estimates of blood velocity and flow rate. We now turn our attention to previously proposed models of
retinal blood flow.

Article submitted to Royal Society

Bernabeu, Jones, Nielsen, et al.
6

Measurement
Arterial diameter
Venous diameter
Mean arterial velocity
Mean venous velocity
Arterial flow rate
Venous flow rate
Arterial flow rate

Venous flow rate

Value
∼ 57 µm
∼ 62 µm
∼ 25 mm s−1
∼ 24 mm s−1
∼ 3.9 µl min−1
∼ 4.8 µl min−1
(0.40 ± 0.04) µl min−1 ,
(0.55 ± 0.06) µl min−1 ,
(0.50 ± 0.05) µl min−1 ,
(0.48 ± 0.05) µl min−1 ,
(0.40 ± 0.04) µl min−1 and
(0.45 ± 0.06) µl min−1
(0.45 ± 0.04) µl min−1 ,
(0.62 ± 0.06) µl min−1 ,
(0.60 ± 0.04) µl min−1 ,
(0.59 ± 0.05) µl min−1 ,
(0.38 ± 0.03) µl min−1 and
(0.63 ± 0.06) µl min−1
(60.4 ± 0.7) µm
(69.3 ± 1.3) µm
(28.3 ± 1.4) mm s−1
(26.3 ± 1.2) mm s−1
up to 28 µm
5 µm to 6 µm
(91.23 ± 11.80) µm
(24.15 ± 1.50) mm s−1
(6.83 ± 1.75) µl min−1
(69.83 ± 3.52) µm
(46.43 ± 1.42) mm s−1
(6.42 ± 0.72) µl min−1

Table 1. Experimentally observed values of vessel diameter and flow rate in different parts of the murine and human adult retina as reported in
various publications. BTBR and C57BL/6J are two common mouse strains used as models of human disease. Values given as mean ± standard
error of the mean, when possible. The values reported by Zhi et al. (2012) are measured at five different arteries/veins and averaged over three
independent measurements. Significant inconsistencies are found across the surveyed literature: a) Wright et al. (2012) and Wang et al. (2007)
measured flow rates approximately one order of magnitude higher than Zhi et al. (2012), b) the arterial diameter measured in Ninomiya and
Inomata (2006) is substantially lower than in Wright et al. (2012).

Species
30-week-old C57BL/6J male mice

11–12-week-old C57BL/6 male mice

Authors
Wright et al. (2012)

Wright et al. (2009)

Four-month-old mice (ex vivo)

22-week-old BTBR female mice

Ninomiya and Inomata (2006)

Adult human

Zhi et al. (2012)

Wang et al. (2007)

Arterial diameter
Venous diameter
Mean arterial velocity
Mean venous velocity
Arterial diameter
Capillary diameter
Arterial diameter
Arterial peak velocity
Arterial flow rate
Venous diameter
Venous peak velocity
Venous flow rate

Article submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

7

Table 2. Retinal peak velocities measured in CD-1 mice by Brown et al. (2005). σ:
standard deviation, N: number of samples.
Velocity (cm s−1 )
0.32
0.71
1.22
0.58
4.62
2.70
5.26
4.30
3.53
2.79
4.35

Age
P0
P1
P2
P3
P4
P5
P6
P7
P8
P10
P12

σ (cm s−1 )
0.09
0.29
0.29
0.07
1.09
0.94
1.79
0.07
1.73
0.14
0.52

N
5
7
6
6
5
5
5
6
7
5
5

Table 3. Values of mean arterial pressure (MAP) and intraocular pressure (IOP)
reported in the literature for different species.
Authors
Wright and Harris (2008)
Hardy et al. (1994)

Species
16-week-old mice
1–3-day-old piglets

MAP (mmHg)
68.2 ± 2.0
70 ± 6

IOP (mmHg)
11.6 ± 0.4
13 ± 2

Table 4. Blood viscosity as function of shear rate in mice.
Authors
Vogel et al. (2003)

Animals
4–7-month-old C57Bl/6 mice

Windberger et al. (2003)

4–8-month-old BALB/c mice

Shear rate (s−1 )
2
5
11
23
45
90
225
450
0.7
2.4
94

Viscosity (mPa s)
18.94 (average, n = 11)
13.33
10.52
8.07
6.31
5.96
4.91
3.85
13.36 (median, n = 37)
10.56
4.87

(c) Previous modelling and simulation studies
In one of the earliest works on retinal haemodynamics modelling and simulation, Ganesan et al. (2010) state that “although a relatively good understanding
of the retinal anatomy and vascular network has been developed through extensive
studies [...] there is a complete lack of numerical modeling of retinal circulation in
the literature”. In the same study, an image-based network model of the retinal vasculature was developed. The location and length of non-capillary vessel segments
was extracted from confocal microscopy images of flat-mounted mouse retinas and
a rule-based network model used to approximate the structure of the capillary bed.
The haemodynamics were greatly simplified by considering vessel segments to be
Article submitted to Royal Society

8

Bernabeu, Jones, Nielsen, et al.

straight with piece-wise constant radius and flow to be laminar (i.e. a 1D network
model), therefore neglecting complex fluid patterns that may appear in curved vessels even at low Reynolds numbers (Siggers and Waters, 2005). Their results show
that wall shear stress in the capillaries stays mainly in the 4 Pa to 11 Pa range with
values as high as 20 Pa. These magnitudes are substantially higher and with a much
greater spread than those reported in the main retinal arteries and veins. Prior to
this work, Liu et al. (2009) also developed an image-based retinal flow model for
the study of oxygen transportation in the retina. In this case, only a subset of
the retinal vasculature (i.e. an artery and a number of branching arterioles) was
reconstructed from a healthy human fundus camera image. The two-dimensional
steady-state Navier-Stokes equations were solved in the domain. The model was
used to predict pressure drops and oxygen saturation distribution.
More recently, Chen et al. (2012) developed a mathematical model of blood flow
in the zebrafish larvae midbrain vasculature (another typical model for the study
of vascular development). The morphology of the vessel network was recovered
from in vivo images at different developmental stages. Haemodynamics were also
modelled using a laminar flow in straight circular pipe simplification. Both steady
and pulsatile flow were compared with little difference in overall dynamics. The
flow model was in turn coupled to a phenomenological model of changes in vessel
diameter as a function of shear stress (without any explicit mechanism of endothelial
cell migration). The resulting coupled model was used to predict vessel pruning in
several zebrafish larvae midbrain vasculature samples with a reported 75% accuracy.
These results support the hypothesis that wall shear stress is a major factor in vessel
pruning during angiogenesis.
Finally, Watson et al. (2012) developed a comprehensive model of murine retinal
angiogenesis including cell migration during sprouting, blood flow, oxygen distribution, and the main chemotactic gradients involved in vessel development and
pruning. In their work, vascular pruning was mainly driven by the downregulation
of growth factors but no explicit mechanobiological mechanisms were considered.
Blood flow simulation was also performed based on the 1D network simplification
described above.
The articles cited in this section demonstrate increasing interest in the modelling
and simulation of retinal haemodynamics. Computational models have been developed for human, mouse, and zebrafish retinal vasculature. A common application
is the study of vascular development dynamics (angiogenesis in particular). In our
opinion, the previous works, although seminal, share one or more of the following
limitations: i) the choice of imaging modality only allows the recovery of a subset
of the retinal vasculature, ii) flow dynamics are greatly simplified by the use of 1D
or 2D approximations, and iii) the computer code developed is, to the best of our
knowledge, not freely available. In the current work, we aim at developing an open
source computational workflow for the generation of high resolution estimates of the
haemodynamic forces experienced by murine retinal vasculature based on confocal
microscope images. The following section describes the methodology employed.

Article submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

9

3. Methods
(a) Image processing and 3D model reconstruction
The preparation of retinal vascular plexus samples for imaging and analysis has been previously described in Franco et al. (2008). Briefly, plexus samples
were collected from five and six-day-old wildtype mouse pups and fixed with 2%
paraformaldehyde in phosphate buffered saline (PBS) for 5 hours at 4◦ C, thereafter retinas were dissected in PBS. Blocking/permeabilisation was performed using Claudio’s blocking buffer (CBB), consisting of 1% FBS (Gibco), 3% BSA
(Sigma), 0.5% triton X100 (Sigma), 0.01% Na deoxycholate (Sigma), 0,02% Na
Azide (Sigma) in PBS pH=7.4 for 2-4 hours at 4◦ C on a rocking platform. Samples were stained with endothelial luminal marker (ICAM2) and incubated at the
desired concentration in 1:1 CBB:PBS at 4◦ C overnight in a rocking platform. Finally, retinas were mounted on slides using Vectashield mounting medium (Vector
Labs, H-1000) and imaged with a Carl Zeiss LSM780 scanning confocal microscope
(Zeiss). Manual preprocessing of the plexus image was performed with Photoshop
CS5 (Adobe) in order to remove major imaging artefacts, allowing a simple thresholding method to be used to produce a binary image of the entire retinal plexus
segment. The binary image was skeletonised using a MATLAB (The MathWorks,
Inc.) interface† to the Voronoi tessellation algorithm implemented in the QHull library (Barber et al., 1996). The radius at each Voronoi vertex was calculated from
the maximum inscribed circle in two-dimensions.
Based on the image skeleton and computed radii, a three-dimensional triangulation of the plexus luminal surface was generated by assuming vessel circular cross
section with the VTK (Schroeder et al., 2003) and VMTK (Antiga et al., 2008)
libraries. This simplification is made based on our own histological analysis and
on the observations of Feke et al. (1989), who cite histological evidence of retinal
arteries being circular in cross-section while retinal veins exhibit a higher tendency
towards flattening. The volume contained within the surface was discretised as a
regular grid in order to generate the computational domain necessary for simulation (see next section for details). All the scripts developed are freely available at
https://github.com/UCL/BernabeuInterface2014.

(b) Simulation setup
Let Ω be a three-dimensional domain with boundary ∂Ω. The CFD package
HemeLB (Mazzeo and Coveney, 2008) was used to solve numerically the Navier
Stokes equations for generalised Newtonian incompressible fluids. For x ∈ Ω̄ and
time t ∈ [0, T ]:


ρ

∇·v =0,

= −∇P + ∇ · T ,

∂v
+ v · ∇v
∂t

(3.1)
(3.2)

† http://www.mathworks.co.uk/matlabcentral/fileexchange/27543-skeletonization-using-voronoi

Article submitted to Royal Society

10

Bernabeu, Jones, Nielsen, et al.

where ρ is the density, v(x, t) is the velocity vector, P (x, t) is the pressure, T(x, t)
is the deviatoric part of the stress tensor
Tij = 2ηSij ,


1 ∂vj
∂vi
Sij =
+
,
2 ∂xi
∂xj

(3.3)
(3.4)

and η is the dynamic viscosity which may depend on the shear rate γ̇, i.e. η(γ̇),
p
(3.5)
γ̇ = 2Sij Sij ,
where i, j = 1, 2, 3 and summation over repeated indices is assumed. Note that in
the case of Newtonian fluids, η(γ̇) = η = const.
Let ∂Ωw,i,o be the wall, inlet, and outlet portions of the domain boundary,
respectively, such that ∂Ω = ∂Ωw ∪ ∂Ωi ∪ ∂Ωo . Equations (3.1)–(3.2) are closed
with the following initial condition
v(x, 0) = 0 ,

x ∈ Ω̄ ,

(3.6)

and boundary conditions
v=0,
η
P n̂ − ∇v · n̂ = Pi n̂ ,
ρ
η
P n̂ − ∇v · n̂ = Po n̂ ,
ρ

x ∈ ∂Ωw ,

(3.7)

x ∈ ∂Ωi ,

(3.8)

x ∈ ∂Ωo ,

(3.9)

(i.e. a pressure drop problem as formulated in Heywood et al. (1996); Formaggia
et al. (2009)) where n̂(x), x ∈ ∂Ω, is the boundary normal vector and Pi,o (t)
are the pressures at the inlet and outlet, respectively. HemeLB uses the latticeBoltzmann (LB) algorithm (see Appendix A for a brief introduction) and runs
efficiently on large scale high performance computing resources (Groen et al., 2013).
HemeLB’s source code is available under LGPL licence and can be downloaded
from http://ccs.chem.ucl.ac.uk/hemelb. Simulations were run either locally
(Section 4a) or on up to 5040 cores of ARCHER, UK National Supercomputing
Service (Section 4c, Appendix B, and Appendix C)
Experimental measurements were used to derive a functional form for η(γ̇) based
on the Carreau-Yasuda (CY) mathematical model (see e.g. (Boyd et al., 2007)):
η(γ̇) = η∞ + (η0 − η∞ ) [1 + (λγ̇)a ]

n−1
a

,

(3.10)

where a, n, and λ are empirically determined to fit a curve between regions of
constant η∞ and η0 . This model defines three rheological regimes: a Newtonian
region of viscosity η0 for low shear rate, followed by a shear-thinning region where
η decreases with γ̇; finally, a second Newtonian region of viscosity η∞ is defined for
high shear rates. Equation (3.10) was fitted to the data in Table 4 with the least
squares algorithm implemented in the gnuplot graphing utility† giving the following
results: η0 = 14.49 mPa s, η∞ = 3.265 mPa s, λ = 0.1839 s, a = 2.707, n = 0.4136.
Figure 2 plots the data fit.
† Gnuplot 4.6.3, http://www.gnuplot.info

Article submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

11

Figure 2. Reported values of murine blood viscosity for different shear rates and
Carreau-Yasuda (CY) model fit.

One of the many challenges when simulating blood flow in open domains — such
as the subset of retinal vasculature that we present in Figure 5 — is the impact of
the choice of inlet/outlet boundary conditions on the simulated haemodynamics.
Ideally, one would use experimental measurements of flow rate and/or pressure in
order to close the system. We could not obtain these data experimentally and relied
on the data surveyed in Section 2b. Pressures at the inlet and outlet were set to the
values measured by Wright and Harris (2008) and presented in Table 3: 68.2 mmHg
and 11.6 mmHg, corresponding to the mean arterial pressure of the central retinal
artery and vein, respectively.
The LB algorithm admits a number of different implementations of the no-slip
boundary condition at the walls (see e.g. Lätt and Chopard (2008); Nash et al.
(2014) for surveys). We choose the method proposed by Bouzidi et al. (2001) based
on previous validation work (Nash et al., 2014). In the current work, we perform
further validation as described in the following section. Finally, we initialise the domain to a uniform density fluid at rest. The implications of this choice are discussed
in Section 4b.

(c) Code verification methodology
In this study, we are interested in using HemeLB to simulate blood flow in
a network of vessels of variable diameter, with differences of up to one order of
magnitude (see Table 1). This is a challenging scenario since we must ensure that
the spatial discretisation is fine enough to capture all the features in the capillaries
and resolve flow accurately, while keeping the problem computationally tractable
due to the large number of fluid sites arising from the discretisation of larger vessels.
Furthermore, we are particularly interested in generating accurate estimates of wall
shear stress in vessels that are typically not aligned with the Cartesian grid, which
can be challenging for regular grid based methods. For example, Stahl et al. (2010)
measured shear stress errors of up to 35% in the vicinity of the wall for nonArticle submitted to Royal Society

12

Bernabeu, Jones, Nielsen, et al.

lattice aligned channel flow with the LB algorithm and the so-called bounce-back
implementation of the no-slip boundary condition (Ziegler, 1993).
In the current section parameters are denoted with a tilde when given in lattice
units and without when given in physical units. The LB timestep ∆t is used as a
conversion factor for time and the voxel size ∆x for space such that e.g. diameter
D = D̃∆x. We also introduce the LB relaxation parameter (which controls the
viscosity in the lattice, for more details refer to Chen and Doolen (1998))
τ̃ =

ν̃
1
,
+
2 c̃2s

(3.11)

where ν̃ is the kinematic viscosity in lattice units
ν̃ =

η ∆t
,
ρ ∆x2

(3.12)

and c̃2s = 31 in the version of LB employed.
In the absence of experimental flow measurements to compare against our computer simulations, we propose setting up a set of benchmark simulations that capture the main flow and domain characteristics and compare the results against
known analytical solutions. We will restrict ourselves to the simulation of steady,
Newtonian, and laminar flow in non-lattice aligned cylinders of diameter D̃ ∈ [3, 30]
and length L̃ = 4D̃. The orientation of the cylinder n̂c is chosen pseudorandomly
from the unit sphere, subject to the constraint that n̂c · êi ≤ 0.9, ∀i. The value is
n̂c = (−0.299, 0.382, 0.874)> .

(3.13)

The laminar flow assumption is based on the Reynolds numbers reported in the
literature for microcirculation (e.g. Re = 0.2, 0.05, and 0.0003 for arterioles, venules,
and capillaries, respectively (Popel and Johnson, 2005)). We choose Re = 1 in our
validation. The steady and Newtonian flow assumptions are made to simplify the
analytical solution of the benchmarks considered. Their implications are discussed
in Section 4a. Finally, a parabolic velocity profile with maximum velocity
ṽmax = −

ν̃Re
n̂c ,
D̃

(3.14)

is imposed at the inlet (Ladd, 1994).
The purpose of our validation study is two-fold. First, to characterise the accuracy of the recovered haemodynamics as a function of the number of fluid sites
across a given vessel. Second, to evaluate the accuracy of the computed wall shear
stress given our choice of implementation of the no-slip boundary condition (Bouzidi
et al., 2001), which has — to the best of our knowledge — not been done before.
In our first experiment, we compare the volumetric flow rate
ZZ
q=
v · n̂ dS
(3.15)
S

integrated over a lattice aligned cross-sectional plane (defined by point (0, 0, 0)>
and plane normal (0, 0, 1)> ) with the analytical solution of Hagen-Poiseuille flow
in an infinite cylinder
|vmax |πD2
q? =
,
(3.16)
8
Article submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

13

for a range of values of diameter D̃ ∈ [3, 30] and relaxation time τ̃ ∈ {0.6, 0.8, 1, 1.2, 1.4}.
In the second experiment, we compare T with the — appropriately rotated —
analytical solution of the Hagen-Poiseuille shear stress tensor T0 in a cylinder of
axis ê3 = (0, 0, 1)> and radius R = D/2 assuming flow in the positive direction
of the cylinder axis. For a given point in the domain x = (x1 , x2 , x3 )> such that
x21 + x22 = r2 , r ∈ [0, R], it can be shown that


0
0
-2|vmax |ηR−2 x1
T0 (x) = 
0
0
-2|vmax |ηR−2 x2  . (3.17)
−2
−2
-2|vmax |ηR x1 -2|vmax |ηR x2
0
The tensor rotation is defined by the matrix R = (r̂1 r̂2 r̂3 ),
r̂1 = n̂c ,
n̂c × ê3
,
r̂2 =
||n̂c × ê3 ||
r̂3 = r̂1 × r̂2 ,

(3.18)
(3.19)
(3.20)

where || · || is the magnitude of a vector, such that
T? = RT0 R> .

(3.21)

4. Results and discussion
(a) Code verification
Figure 3 plots the relative error in the simulated flow rate
q =

q? − q
q?

(4.1)

as a function of the cylinder diameter. We observe how the choice of τ̃ greatly affects
the accuracy of the simulated haemodynamics. In agreement with similar analyses
in the literature (see e.g. Krüger et al. (2009)), the error is larger for values close to
the stability threshold of τ̃ = 0.5 and for values greater than 1. A region of excellent
accuracy is located around τ̃ = 0.8. In that case, the relative error q stays below
3% for all the values of D̃ studied. The ability to simulate correct flow dynamics
in channels with only a few lattice sites across has gained the LB algorithm wide
acceptance for the simulation of flow in complex domains and porous media (Aidun
and Clausen, 2010).
Figure 4 plots, for a range of values of D̃, the computed and analytical solutions
of the shear stress tensor (T and T? ) as well as the associated relative error
T =

||T? − T||F
,
||T? ||F

where || · ||F is the Frobenius norm of an m × n matrix A = [aij ]
v
uX
n
um X
||A||F = t
|aij |2 .
i=1 j=1

Article submitted to Royal Society

(4.2)

(4.3)

14

Bernabeu, Jones, Nielsen, et al.

0.12

τ̃
τ̃
τ̃
τ̃

0.10

q

0.08

= 0.6
= 0.8
=1
= 1.2

0.06
0.04
0.02
0.00
2

4

6

8

10

12

14

16

18

20

D̃
Figure 3. Hagen-Poiseuille flow in an inclined cylinder. Relative error on the computed
flow rate as a function of vessel diameter D̃ and lattice-Boltzmann (LB) relaxation time
τ̃ . For τ̃ = 0.8 the total error is kept below 3% even for cylinders with just 3 lattice sites
across. These results confirm the suitability of the LB algorithm for the simulation of
flow in sparse geometries and porous media. The lines are a guide to the eye and bear no
physical meaning.

The choice of error norm in (4.2) ensures that error contributions given by individual
components cannot compensate, hence being sensitive to rotation errors. We are
mainly interested in the accuracy of the shear stress calculation in the vicinity of
the vessel wall. Therefore, the results are presented for the subdomain defined by
all the lattice sites with r ∈ [0.8R, R].
It can be observed how the recovered shear stress follows the expected pattern
of monotonic increase from zero at the cylinder axis (results not shown here) to
its maximum value at the wall. However, a certain deviation exists compared to
the analytical solution. For values of D̃ ≥ 7, we observe a clear pattern of shear
stress being overestimated in the [0.8R, 0.9R] region while being underestimated
in [0.9R, R]. Panels 4(b), 4(d), and 4(f) quantify this error. We observe how the
largest error in the domain always occurs at the cylinder wall and that it decreases
as D̃ increases: from around 17% for D̃ = 5 to around 7% for D̃ = 15. These results,
with the Bouzidi et al. (2001) implementation of the no-slip boundary condition,
represent a substantial improvement over the 35% error reported by Stahl et al.
(2010) with the bounce-back method and D̃ = 20, confirming the superiority of the
former algorithm. More importantly, the wall shear stress is consistently underestimated and the error decreases for D̃ ≥ 7. Therefore, one could implement an a
posteriori correction of the computed shear stress tensor based on this knowledge.
Such a development is beyond the scope of the current paper and will be developed
as part of a future study.
It is interesting to note that the errors in flow rate and shear stress are, to a large
extent, decoupled. This is due to the characteristics of the LB algorithm as detailed
Article submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

0.00144
0.00140

0.18

||T||F

0.16

||T?||F

0.14

0.00132

0.12

0.00128

0.10

T

Pa

0.00136

0.00124

0.08
0.06

0.00120

0.04

0.00116

0.02

0.00112
0.8R

0.9R
r

R

0.8R

(a) D̃ = 5, analytical vs. computed solution.
0.00144
0.00140

0.16

||T?||F

0.14
0.12

0.00128

0.10

0.00124

0.08
0.06

0.00120

0.04

0.00116

0.02

0.00112
0.8R

0.9R
r

R

0.8R

(c) D̃ = 7, analytical vs. computed solution

0.00140

R

0.18

||T||F

0.16

||T?||F

0.14

0.00132

0.12

0.00128

0.10

0.00124

0.08
0.06

0.00120

0.04

0.00116
0.00112
0.8R

0.9R
r

(d) D̃ = 7, relative error.

T

Pa

0.00136

R

0.18

||T||F

0.00132

0.00144

0.9R
r

(b) D̃ = 5, relative error.

T

Pa

0.00136

15

0.02
0.9R
r

R

(e) D̃ = 15, analytical vs. computed solution.

0.8R

0.9R
r

R

(f) D̃ = 15, relative error.

Figure 4. Hagen-Poiseuille shear stress in inclined cylinders of D̃ ∈ {5, 7, 15}. Norm of
the analytical and computed stress tensors (panels 4(a), 4(c), 4(e)) and relative error
between them (panels 4(b), 4(d), 4(f)). Results are presented for every lattice site with
radius r ∈ [0.8R, R]. Agreement between computed and analytical solution improves with
increasing D̃ and it is always better in the vicinity of the cylinder wall. These results, with
the Bouzidi et al. (2001) implementation of the no-slip boundary condition, represent
a substantial improvement over the 35% error reported by Stahl et al. (2010) with the
bounce-back method and D̃ = 20.

in Appendix A. In our work we will define two rules for the accurate simulation
of blood flow in our network of vessels of variable diameter. First, a minimum
diameter of D̃ = 3 will be enforced throughout the domain. This will ensure that,
for τ̃ = 0.8, the general flow patterns produced are accurate. Second, values of
shear stress in regions of interest will only be considered valid if D̃ ≥ 7. The error
Article submitted to Royal Society

16

Bernabeu, Jones, Nielsen, et al.

estimates reported in Figure 3 will be taken into account in the presentation of our
results.
Finally, we turn our attention to the implications of the use of a generalised
Newtonian rheology model in our retinal flow simulations. In the current section
we observed that τ̃ = 0.8 minimises the error in the flow rate recovered. In the
case of Newtonian fluids, regardless of the value of viscosity being simulated, one
∆t
can choose ∆x
2 such that Equation (3.11) yields the desired value of τ̃ (this will
obviously have an impact on computational cost). For generalised Newtonian fluids,
τ̃ becomes a function of γ̇ and will take values in the range [τ̃∞ , τ̃0 ]. In the rheology
0
= 4.92, which would lead to impractical values
model presented in Section 3b, νν∞
∆t
of τ̃∞ unless τ̃0 is chosen small enough. We will therefore choose the coefficient ∆x
2
such that τ̃∞ = 0.6 and τ̃0 = 0.992. Based on the results in Figure 3, this choice
will yield an error of less than 4% in the flow rate recovered for D̃ ≥ 7.
(b) Model reconstruction
The vascular plexus of wildtype retinas was stained with the luminal membrane
marker ICAM2 and images were acquired using a confocal microscope as described
in Section 3a. Figure 5(a) shows a region of interest in one of the imaged retinas.
It contains, on either side, two arterial segments coming from the optic disc and
connecting with a segment of a retinal vein (centre of the image) through a dense
capillary network. It can be appreciated how the network is more mature (e.g. vessel
identity, branching patterns) closer to the optic disc (bottom of the image) while
its structure is much more primitive and less remodelled in the periphery closer to
the sprouting front (top of the image). Figure 5(b) presents the results of the image
segmentation process. The algorithm described in Section 3a is used to first create a
binary mask separating the luminal area and background tissue and second extract
the network skeleton and radii. The latter are used to reconstruct the 3D luminal
surface under the assumption of vessel circular cross-section (see Section 3a for a
discussion). Figure 5(c) shows the reconstructed surface. We refer to this model as
P6A. Figure 6 presents luminal surface binary masks for three additional P5 and
P6 retinal plexuses. The same reconstruction algorithm is applied and the resulting
models are referred to as P5A, P5B, and P6B, respectively.
Figure 7 plots a network diameter histogram (in terms of total distance covered by vessel segments of a given diameter) for models P5B and P6B. The largest
diameter in the network are Dmax = 34 µm and 40 µm, respectively, which occur
along the retinal vein. The artery segments have diameters of up to 16 µm, with
larger diameters closer to the optic disc. The bulk of the capillary bed has diameters approximately in the range 2 µm to 10 µm, with a reduced amount of vessels
with smaller diameter. These results are substantially lower than the in vivo measurements presented in Table 1. In addition, a small number of capillaries have
diameters approaching 0 µm. To some extent, these discrepancies are expected as
we are measuring the diameter of the internal luminal surface with extreme precision (unlike some of the works cited where only a generic measure of vessel calibre is
given), including in our measurements vessel segments that appear to be undergoing
regression (hence in the process of closing up). The diameter measured for the main
arteries are, however, in better agreement with the ex vivo measurements obtained
from corrosion casts by Ninomiya and Inomata (2006). Therefore, we cannot exArticle submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

(a)

17

(b)

(c)

Figure 5. Subset of a wildtype P6 retinal plexus used to reconstruct one of our retinal
blood flow models, namely P6A model. The original microscope image is segmented and
the network skeleton and segment radii are computed. Based on these values, a 3D volume
is reconstructed assuming vessels of piece-wise constant radius. (a) Original image. (b)
Segmented image. (c) Reconstructed surface.

clude that sample preparation and fixation protocols contribute to vessel shrinkage.
Finally, we fit a lognormal probability distribution function to each histogram and
use the distribution mode as an estimate of the typical capillary diameter (under
the assumption that capillaries are the most common vessel type in the network).
We observe a reduction in the typical capillary diameter between day five (5.51 µm
and 5.78 µm) and six (4.44 µm and 5.29 µm). More experiments are required in
order to assess the statistical relevance of these results but the implications of a
systematic decrease in vessel diameter over time are important given that, for a
constant flow rate, wall shear stress is inversely proportional to the third power of
the vessel radius. We plan to explore the relationship between changes in geometry
and haemodynamics as part of a future study.
In order to ensure that 95% of the reconstructed network has D̃ ≥ 7, we choose
the voxel sizes in Table 5 for the discretisation of each model. Appendix B presents
a grid refinement study aimed at confirming that the choice of voxel size leads to
spatially converged solutions.
Article submitted to Royal Society

18

Bernabeu, Jones, Nielsen, et al.

(a) P5A model.

(b) P5B model.

(c) P6B model.

Figure 6. Binary masks defining the luminal surface of three retinal plexuses obtained
at two different stages of development. All plexuses are presented with the area closer
to the optic disc at the bottom of the image and the sprouting front at the top. In all
samples studied, arteries tend to be thinner and have less daughter vessels than veins.
Vessels close to the sprouting front tend to have less well-defined identity with luminal
diameters comparable to arteries/veins. This is particularly noticeable in the P5 samples.
Vessel density is also higher close to the sprouting front in P5 retinas.

(c) Simulations
Due to its kinetic nature, the LB algorithm applied to steady flow problems in
an initially quiescent domain requires the system to be advanced in time in order
to overcome an initial transient. In order to monitor convergence, we evaluate the
following convergence criterion at the end of each timestep t:
maxr ||v(r)t − v(r)t−1 ||
< tol = 10−6 ,
vref

(4.4)

where v(r)t := v(r, t∆t) and vref is a velocity reference value chosen based on the
data summarised in Table 1, i.e. vref = 50 mm s−1 . Only when this condition holds
do we consider the simulation to have reached steady state. More efficient methods
for LB initialisation have been proposed (see e.g. Mei et al. (2006); Caiazzo (2005))
but we will not consider them in the current work since our approach remains
computationally tractable.
Figure 8 presents results of a simulation with the P5B flow model and the
inlet/outlet boundary conditions and rheological properties surveyed in Section 2b.
Velocity magnitude is plotted at the intersection of the model and the z = 0 plane
(Figure 8(a)). Our results show how velocities are larger in the central artery (see
label A), in particular close to the optic disk. We also note that, as the artery
progresses towards the sprouting front, the velocity magnitude decreases rapidly.
Furthermore, it stops being a preferential flow path at the point where it meets
areas of less well established vessel identity close to the sprouting front (see e.g.
B regions). Areas with undefined vessel identity are correlated with homogeneous
velocity distributions (see e.g. C). There exists evidence of a considerable number of
vessels having recently regressed along the path of the artery (see e.g. D branches)
and the more developed first order branches (see e.g. E). The two veins (top and
bottom of the images) present fewer regressing profiles.
Article submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

19

(a) P5B model.

(b) P6B model.

Figure 7. Network diameter histogram showing the aggregated total distance covered by
vessels of a given diameter. Vertical lines indicate the mode of a lognormal probability
distribution fit of each dataset. The values for the models not shown here are 5.51 µm
(P5A) and 5.29 µm (P6A). We use these values as an estimate of the typical capillary
diameter (the most common type of vessel in the network). Capillaries with diameter
approaching 0 µm appear to be undergoing remodelling. Arterial and venular segments
present higher diameters ranging up to 34 and 40 µm, respectively.
Table 5. Voxel sizes employed in the discretisation of the different flow models used in
this work.
Model
P5A
P5B
P6A
P6B
∆x
0.5166 µm 0.5666 µm 0.5 µm 0.4166 µm

Article submitted to Royal Society

20

Bernabeu, Jones, Nielsen, et al.

For a plane with normal n̂, we define the traction vector
t = Tn̂ ,

(4.5)

i.e. the force per unit area acting on that plane. Figure 8(b) plots traction magnitude
||t|| on the model surface (often referred to as wall shear stress magnitude, WSS).
We observe that areas of preferential flow correlate well with the areas experiencing
larger WSS. In contrast, vessels in the sprouting front are under lower magnitudes
of WSS. The model predicts values of WSS larger than 20 Pa, which can be deemed
unphysiological based on the microvasculature WSS measurements reported in the
literature: 14 Pa (Popel and Johnson, 2005), ∼ 20 Pa (Ganesan et al., 2010)) or
∼ 13 Pa computed from the values reported by Wright et al. (2012) (under the
assumption of Poiseuille flow). We believe that the WSS overestimation (mainly
occurring at the central artery and some first order branches) is due to the vessel
shrinkage discussed earlier or other modelling errors.
Figure 9 presents results of a simulation with the P6A flow model and the inlet/outlet boundary conditions and rheological properties surveyed in Section 2b.
Figure 9(a) plots velocity magnitude on the intersection of the model and the z = 0
plane. First of all, it can be appreciated how velocities are larger on arteries, veins,
and first-order vessels branching out from them. Highest peak velocities are around
42 mm s−1 (corresponding to mean velocities of 21 mm s−1 under Poiseuille flow
assumption) and are in good agreement with the measurements by Wright and
colleagues presented in Table 1. Velocity distribution along a given vessel diameter displays the expected parabolic profile with zero velocity at the walls. Areas
in more advanced state of pruning (typically closer to the optic disc, see region
B) tend to present larger velocity magnitudes due to a reduction in vessel density.
An exception to this trend is region C. In this case, we observe a region of very
low flow (similar to the regions found in the less mature vascular plexus towards
the periphery) in an area where pruning should be in a fairly advanced stage. Two
explanations are possible: i) that a recent vessel regression event has drastically
reduced the total flow arriving to the area which in turn will trigger further vessel
regression (similar to what can be observed in region A) or ii) that a vessel segment
connecting the area with the nearby artery was accidentally removed when preparing the sample. In contrast to the optic disc region, areas in the vicinity of the
sprouting front experience lower velocity magnitudes (see region E). Nevertheless,
even in this region, we can already appreciate segments of predominant flow (see
e.g. F regions) rather than a totally homogeneous flow distribution. Taken together
with observations from Chen et al. (2012), we predict that these high-flow vessel
segments are likely to survive the pruning process.
Figure 9(b) plots WSS on the model surface. We observe that areas of preferential flow correlate well with the areas experiencing larger WSS. However, in this
case we do not see a decrease in WSS with increasing vessel order. WSS peaks are
distributed throughout the domain in agreement with the observations by Ganesan
et al. (2010). We also observe a complex distribution of WSS along individual vessel
segments, with changes following local variation in vessel diameter.
Figure 10 plots t on the surface of a subset of the domain (marked with a circle
in Figure 9(b)). Given the redundancy of a loop-like structure of this type and
the distribution of diameters present, it can be assumed that the upper half of the
loop is undergoing regression. This fact is in good agreement with the distribution
Article submitted to Royal Society

Sprouting front

Computational analysis of haemodynamics during angiogenesis




B eD
h
A
De

E

21

h


'$
B

C&%

Sprouting front

(a)

(b)

Figure 8. P5B simulation results: (a) Velocity magnitude plotted on a cross-section along
the z = 0 plane. Velocity shows the expected parabolic profile across the vessel diameter.
Velocity is higher in the artery located in the centre of the domain, in particular close to
the optic disk. Velocity magnitude quickly decreases as the artery progresses towards the
sprouting front and it stops being a preferential flow path at the points where its identity
stops being clearly defined. (b) Wall shear stress (WSS) magnitude plotted on the model
surface. Areas of preferential flow tend to experience highest WSS magnitudes. WSS is
generally low across the domain except for the arterial segment close to the optic disk and
some first order branches. WSS values higher than 20 Pa are considered unphysiological
and the regions experiencing them are coloured in black. Black circles indicate regions of
interest referenced in the manuscript.

of WSS, of much larger magnitude on the bottom section of the loop and vessel
segments upstream and downstream from it.
In summary, the results presented in this section support the idea that vessel
segments undergoing pruning tend to occur in regions of low flow (and hence low
Article submitted to Royal Society

22

Bernabeu, Jones, Nielsen, et al.

Sprouting front

h
F

hF

hA


E





D








C

B

(a)

Sprouting front



(b)

Figure 9. P6A simulation results: (a) Velocity magnitude plotted on a cross-section along
the z = 0 plane. Velocity is higher in arteries, veins, and segments directly branching
from them close to the optic disc. Velocity magnitude is smaller in the sprouting front.
However, vessels of preferential flow already exist in the sprouting front; potentially an
early indicator of which vessels will survive the pruning process. (b) Wall shear stress
(WSS) magnitude plotted on the model surface. Areas of preferential flow tend to experience highest WSS magnitudes. WSS peaks are widely spread across the network. WSS
magnitude tends to be lower at the junctions and many vessel segments present a high-low pattern due to local changes in vessel diameter. WSS values higher than 20 Pa are
considered unphysiological and the regions experiencing them are coloured in black. Black
circles indicate regions of interest referenced in the manuscript.

shear stress). We hypothesise that this process gradually reduces network density
and as a consequence flow increases in the surviving vessel segments. This in turn
prevents further pruning and contributes to vessel maturation.
Article submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

23

Figure 10. Traction vectors (of constant length and coloured according to magnitude) on
the luminal surface of the region of interest highlighted in Figure 9(b). The loop branch
undergoing regression (upper branch) experiences a much lower traction magnitude.

(d ) Limitations of the study
The main limitations of the study are as follows. First, blood was modelled as
a homogeneous fluid rather than a particle suspension. This is likely to have an
impact on the wall shear stress computed in small calibre capillaries. Xiong and
Zhang (2010) studied the changes in haemodynamics induced by the presence of
red blood cells flowing in a simplified model of a microvessel and found up to a
20% increase in the shear stress experienced by the luminal wall. Second, although
blood was modelled as a shear-thinning fluid, other rheological properties such
as the Fåhræus-Lindqvist effect (see Section 2b) were not accounted for. Third,
vessel cross-section was assumed to be circular throughout the domain due to the
lack of spatial information in the z axis. As previously mentioned, there exists
experimental evidence supporting this assumption in retinal arteries but not in
veins (Feke et al., 1989). This will have an impact on the haemodynamics recovered.
Also, despite all our efforts when processing retina samples, we cannot be fully
certain that no distortions in the vascular plexus were introduced. Next, due to the
difficulty of measuring pressure or flow profiles at the model inlets/outlets in vivo
and the absence of suitable data in the literature, only steady state simulations
were performed. We expect flow to be nearly in phase with pressure given the
typical values of Womersley number (defined as the ratio between oscillatory inertial
forces and viscous forces) encountered in retinal circulation (∼ 0.1 according to
Liu et al. (2009)). This makes us confident that flow has time to fully develop
in each cardiac cycle and hence will be well approximated by an instantaneous
pressure gradient. Nevertheless, there will still be substantial variations in wall
shear stress within any given cardiac cycle. Furthermore, the values of mean arterial
pressure (MAP) and intraocular pressure (IOP) used as inlet and outlet boundary
conditions were obtained from adult animals. In the supplementary material of
this paper (Appendix C), we perform a sensitivity analysis of these parameters.
Finally, another source of variation in the predicted haemodynamics are the active
Article submitted to Royal Society

24

Bernabeu, Jones, Nielsen, et al.

and passive mechanical properties of retinal vessels. At the analysed stage, retinal
arteries are already covered with a smooth muscle layer, which might contract/relax
to control local flow (i.e. autoregulation) and therefore have an impact in flow
patterns in downstream vessels.

5. Conclusions
In the current work, we have presented a software pipeline for the creation of
computational blood flow models based on confocal microscope images of the microvasculature. The pipeline has been applied to the development of flow models of
the neonatal mouse retinal vasculature (a common animal model for the study of
vascular development). The different software components used are released under
open source licences.
Using simplified benchmark problems, we have demonstrated the suitability of
the lattice-Boltzmann (LB) algorithm for the simulation of blood flow in sparse and
highly complex vascular networks. Our results indicate that a careful choice of the
LB configuration parameters leads to accurate flow estimates in channels as narrow
as three lattice sites across. Furthermore, we also showed that the implementation
of the no-slip boundary condition proposed by Bouzidi et al. (2001) produces acceptable estimates of wall shear stress. We measured errors of ∼ 10% and ∼ 7%
in channels 7 and 15 lattice sites wide, respectively. Being able to recover correct
haemodynamics even at moderately coarse discretisations is fundamental to keep
the problems under study computationally tractable.
In the study reported here, we investigated changes in haemodynamics during
vascular remodelling. Blood flow models were generated from samples of retinal
plexuses obtained at postnatal day (P) 5 and 6. Our simulations show that, in
both cases, velocity and wall shear stress (WSS) are higher in arteries, veins, and
first-order capillaries closer to the optic disk. However, important differences in the
distribution of velocity and WSS across the domain are observed when comparing
both days (see e.g. Figures 8(b) and 9(b)). On the one hand, P5 simulations show
a very homogeneous distribution of velocity and WSS across the capillary network
with moderately high values only in the vicinity of the optic disk. On the other
hand, simulations with the P6 flow model show a consistently higher and much
more spatially complex distribution of velocity and WSS. Higher values are primarily located in regions in a more advanced state of remodelling (note e.g. the
number of disconnected vessels undergoing regression). In the P6 case, branches of
predominant flow can be also identified in the sprouting front.
We also analysed WSS in segments undergoing regression (see e.g. Figure 10)
and observed vessel pruning occurring in regions of low shear stress. This process
gradually reduces network density (through the removal of redundant segments)
and is likely to lead to an increase in flow in the surviving vessel segments. We
hypothesise that this will contribute to vessel maturation. Our results support the
previously proposed modulation effect that haemodynamic forces have on developmental vascular remodelling (Chen et al., 2012).
The geometrical analysis of the vascular plexuses leads to two possible explanations for the increase in velocity and WSS observed between the P5 and P6 models.
First, the increase may be a direct consequence of the observed decrease in typical capillary diameter (given the inverse relationship between WSS and the third
Article submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

25

power of the vessel radius and assuming that the total flow rate in the retina remains
constant). Second, the progressive reduction in capillary bed density due to vessel
regression may lead to an increase in flow (and hence WSS) in neighbouring vessels. We believe that both effects may play complementary roles in order to create
the WSS gradients hypothesised to be behind vessel regression (Chen et al., 2012).
Further experiments are required in order to determine the relative importance of
each of these effects and fully understand how they interact.
We are currently working on extending the modelling framework to include
tissue mechanics and agent-based cellular modelling. Our goal is to develop an
integrated computational framework for vascular mechanobiology research. In particular, we are interested in modelling the interplay between cellular molecular
regulation and haemodynamic forces during vascular remodelling. Finally, the developed methodology should be applicable to other research domains where small
vascular networks can be imaged but where experimental flow measurements are
difficult to obtain.

Acknowledgements
This work was supported by: Cancer Research UK; the Lister Institute of Preventive Medicine; the Leducq Transatlantic Network ARTEMIS; the UK-Israel Initiative BIRAX; EPSRC grants “2020 Science” (http://www.2020science.net/,
EP/I017909/1), “Large Scale Lattice Boltzmann for Biocolloidal Systems” (EP/I034602/1),
and “UK Consortium on Mesoscale Engineering Sciences (UKCOMES)” (EP/L00030X/1);
and the EC-FP7 projects CRESTA (http://www.cresta-project.eu/, grant no.
287703) and MAPPER (http://www.mapper-project.eu/, grant no. 261507). CAF
is supported by a Marie Curie Post-doctoral Fellowship of the European Commission FP7 People framework. This work made use of the HECToR and ARCHER
UK National Supercomputing Services (http://www.archer.ac.uk) (under EPSRC
grants EP/I017909/1 and EP/L00030X/1). The authors acknowledge the support
of the UCL Research Software Development Service (RSD@UCL) in the completion
of this work.

Appendix A. An introduction to the lattice-Boltzmann
algorithm
In the current section we present a brief introduction to the lattice-Boltzmann (LB)
algorithm. The interested reader can refer to well-cited references such as Chen and
Doolen (1998); Aidun and Clausen (2010) for a more in-depth presentation and
analysis. LB operates at a mesoscopic level, simulating the evolution of a discretevelocity approximation to the one-particle velocity distribution functions of the
Boltzmann equation of kinetic theory, {fi (r, t)}. Computations are performed on
a regular lattice discretisation of Ω̄w , with grid spacing ∆x. The set of velocities
{ci } is chosen such that the distances travelled in one timestep (∆t), ei = ci ∆t,
are lattice vectors. When one only wishes to reproduce Navier-Stokes dynamics, the
set is typically a subset of the Moore neighbourhood, including the rest vector. For
3D simulations, the most commonly used sets have 15, 19 and 27 members. In the
current work, we employ the three-dimensional 15 velocity LB lattice (D3Q15).
Article submitted to Royal Society

26

Bernabeu, Jones, Nielsen, et al.

Evolving the distribution functions in time involves two main steps. The first is
known as the collision step, which relaxes the distributions towards a local equilibrium (the post-collisional distributions is often denoted as fi? ):
fi? (r, t) = fi (r, t) + Ω̂(fi (r, t)) ,

(A 1)

where Ω̂ is a the collision operator. The second is known as the streaming step,
where {fi? } are propagated along the lattice vectors to new locations in the lattice,
defining the distribution functions at the next timestep:
fi (r + ci ∆t, t + ∆t) = fi? (r, t) .

(A 2)

In the current work, we employ the lattice Bhatnagar-Gross-Krook (LBGK)
collision operator, which approximates the collision step as a relaxation process
towards a local equilibrium,
Ω̂(fi ) = −

(fi − fieq )
∆t ,
τ

(A 3)

where τ is the relaxation time. This can be shown, through a Chapman-Enskog
expansion (see, e.g., Ladd and Verberg (2001); Chen and Doolen (1998)), to reproduce the Navier-Stokes equations in the quasi-incompressible limit with errors
proportional to the lattice Mach number squared. The kinematic viscosity ν is given
by
ν = c2s (τ − ∆t/2) ,
(A 4)
where
∆x
,
cs = √
3∆t

(A 5)

is the speed of sound in the D3Q15 lattice. Equation (A 3) can be extended to
simulate generalised Newtonian flows by locally varying τ according to (A 4) and
an empirical characterisation of the fluid viscosity (e.g. Equation (3.10)).
For the equilibrium distribution, we use a second-order (in velocity space) approximation to a Maxwellian distribution


ṽ · ṽ
ci · ṽ (ci · ṽ)2
fieq (ρ̃, ṽ) = ρ̃wi 1 +
+
−
,
(A 6)
c˜s 2
2c˜s 4
2c˜s 2
where the weights wi and speed of sound cs depend on the choice of velocity set.
Other choices of Ω̂ exist. Finally, the interested reader can refer to previous work
by the authors (Nash et al. (2014) and references therein) on the implementation
of wall and open boundary conditions.
The macroscopic density ρ(r, t) and velocity v(r, t) at a fluid site can be calculated from the distribution functions by
X
ρ̃ =
fi ,
(A 7)
i

ρ̃ṽ =

X
i

Article submitted to Royal Society

fi ci .

(A 8)

Computational analysis of haemodynamics during angiogenesis

27

The macroscopic pressure is related to the density by the ideal gas law
P̃ = ρ̃c̃2s .

(A 9)

A notable characteristic of the LB algorithm is that the shear rate tensor S can be
computed at any lattice site from local information only. Let fineq = fi − fieq be
the non-equilibrium part of the distribution function. It can be shown that (see e.g.
Krüger et al. (2009)):
S̃ ≈

−1 X neq
f ci ci .
2τ̃ c̃2s ρ̃ i i

(A 10)

The main advantage of this approach is that S (and therefore γ̇, T, and t) can
be evaluated locally without the need of accessing simulation results generated at
neighbouring lattice sites in order to approximate spatial derivatives of the velocity
field. This makes parallel performance independent of the choice of rheology model
or the need of computing wall shear stress (a variable of primary importance in our
work).

Appendix B. Voxel size convergence analysis
In this section we perform a grid refinement study of the simulations presented
in Section 4c. We want to assess whether the characterisation of the discretisation
error done in Section 4a with a simplified domain and rheology model remains valid
with more complex vessel networks and rheology models.
In order to confirm that the choice of ∆x in Table 5 leads to spatially converged
solutions, we generated a set of increasingly finer discretisations of the P6A flow
model. The finest discretisation that remained computationally tractable for us had
∆xr = 0.1875 µm. We used the results of this simulation as a reference solution and
compared it against the results obtained with ∆xi ∈ {1.0, 0.5, 0.25} µm. Note that
∆t had to be modified according to (3.11)–(3.12) in order to keep τ̃∞,0 constant
across all simulations (see Section 4a for more details). We define the error in the
velocity field for a given discretisation ∆xi as
i

i

i

r

r

vi (r∆x ) = v∆x (r∆x , tconv ) − v∆x (r∆x , tconv )

(B 1)

where tconv corresponds, in each of the simulations, to the time where steady flow
convergence has been achieved according to the stop criteria in Equation 4.4 and
r
i
r∆x is the lattice site in the reference discretisation closest to r∆x . Finally, we
use the root mean square of the error scaled by the predicted velocity range as our
measure of error
rX
vi · vi
2ui = √

r

N maxr ||v∆xr ||

(B 2)

where N is the number of lattice sites in the ∆xi discretisation or a subset of it.
Figure 11 plots the values of 2ui computed on the subset of the P6A flow model
presented in Figure 9(a) (i.e. lattice sites on the z = 0 plane) for the choice of
Article submitted to Royal Society

28

Bernabeu, Jones, Nielsen, et al.

Figure 11. Velocity error residual 2ui on simulations with the P6A flow model discretised
with ∆xi = 1.0 µm, 0.5 µm and 0.25 µm (red markers). In order to keep the analysis
computationally tractable, 2ui is computed with results obtained at the lattice sites located
on the z = 0 plane only (this is the same subset of the results presented in Figure 9).
The lines are guides to the eye showing first-order (dashed) and second-order (solid)
convergence.

∆x discussed above. The results show the LB configuration employed (i.e. lattice,
collision operator, and boundary condition implementation) displays second-order
convergence behaviour. WSS (results not presented here) shows similar convergence
trends. These findings are in agreement with results previously published by the
authors (Nash et al., 2014). Finally, our grid refinement study shows how the choice
of voxel sizes in Table 5 leads to sufficiently spatially converged solutions. In the
case of the P6A flow model, the results generated with ∆x = 0.5 µm have a relative
error of only 2ui ≈ 0.005 when compared with the results obtained with finest
discretisation of the model that remained computationally tractable.

Appendix C. Inlet/outlet configuration sensitivity analysis
In this section we investigate the robustness of the results in Section 4c with regards
to the choice of inlet and outlet boundary conditions. This is motivated by the fact
that the choice of MAP and IOP in our flow model is based on data obtained from
the literature (see Table 3) rather than directly measured from the animals used
in the study. First, we want to investigate whether the general patterns of flow
are affected by moderate changes in OPP (defined as the difference between MAP
and IOP). Second, we want to quantify the relationship between OPP and blood
velocity in the domain.
Figure 12 repeats the visualisation in Figure 9(a) for simulations with OPP
10 mmHg smaller and larger. The colour scale has been adjusted to range from
1 mm s−1 to the largest velocity in the domain in each case. It can be appreciated
Article submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

29

(a) 45 mmHg

(b) 65 mmHg

Figure 12. Simulation results: velocity magnitude plotted on a cross-section along the
z = 0 plane for OPP values of 45 mmHg and 65 mmHg. The logarithmic colour scales have
been adjusted to range from 1 mm s−1 to the largest velocity in the domain. Branches of
predominant flow and velocity gradients remain fairly constant despite moderate changes
in OPP when compared to Figure 9(a).

how the general patterns of flow and velocity gradients are greatly preserved from
those in Figure 9. As expected, the absolute values differ.
Table 6 presents the largest velocity magnitude recovered for a wider range of
OPP values. It can be seen how peak velocity increases by approximately 7.6 mm s−1
for every 10 mmHg increase in OPP. The linear relationship between pressure difference driving flow in the domain and peak velocity indicates that flow occurs in the
Stokes regime (i.e. inertial forces are small compared to viscous forces). This result
is expected given the typical Reynolds numbers reported in the literature for microcirculation (e.g. Re = 0.2, 0.05, and 0.0003 for arterioles, venules, and capillaries,
respectively (Popel and Johnson, 2005)).
We conclude that the inlet/outlet boundary conditions play a secondary role in
determining the flow patterns in the network, hence no modulation in the WSS experienced by the endothelium is expected via changes in OPP (this does not include
changes due to autoregulation of the vessels themselves). Changes in vessel network
geometry are expected to be the main drivers behind haemodynamic reorganisation
Article submitted to Royal Society

30

Bernabeu, Jones, Nielsen, et al.

Table 6. Peak velocity in the domain as a function of the ocular perfusion pressure
(OPP).
OPP (mmHg)
25
35
45
55
65

Domain peak velocity (mm s−1 )
18.97
26.65
34.31
41.96
49.61

during development. This observation supports the idea that regression of a given
vessel segment will affect the flow patterns in nearby segments, leading to potential
increases in WSS that may contribute to vessel maturation.

References
Aidun, C. K. and Clausen, J. R. (2010). Lattice-Boltzmann Method for Complex
Flows. Annual Review of Fluid Mechanics, 42(1):439–472.
Alm, A. and Bill, A. (1973). Ocular and optic nerve blood flow at normal and
increased intraocular pressures in monkeys (Macaca irus): a study with radioactively labelled microspheres including flow determinations in brain and some
other tissues. Experimental Eye Research, 15(1):15–29.
Antiga, L., Piccinelli, M., Botti, L., Ene-Iordache, B., Remuzzi, A., and Steinman, D. A. (2008). An image-based modeling framework for patient-specific
computational hemodynamics. Medical & biological engineering & computing,
46(11):1097–112.
Barber, C. B., Dobkin, D. P., and Huhdanpaa, H. (1996). The quickhull algorithm
for convex hulls. ACM Transactions on Mathematical Software, 22(4):469–483.
Bentley, K., Jones, M., and Cruys, B. (2013). Predicting the future: Towards symbiotic computational and experimental angiogenesis research. Experimental cell
research, 319(9):1240–6.
Bernabeu, M. O., Nash, R. W., Groen, D., Carver, H. B., Hetherington, J., Krüger,
T., and Coveney, P. V. (2013). Impact of blood rheology on wall shear stress in
a model of the middle cerebral artery. Interface Focus, 3(2):20120094.
Bouzidi, M., Firdaouss, M., and Lallemand, P. (2001). Momentum transfer of a
Boltzmann-lattice fluid with boundaries. Physics of Fluids, 13(11):3452.
Boyd, J., Buick, J. M., and Green, S. (2007). Analysis of the Casson and CarreauYasuda non-Newtonian blood models in steady and oscillatory flows using the
lattice Boltzmann method. Physics of Fluids, 19(9):093103.
Brown, A. S., Leamen, L., Cucevic, V., and Foster, F. S. (2005). Quantitation of
hemodynamic function during developmental vascular regression in the mouse
eye. Investigative ophthalmology & visual science, 46(7):2231–7.
Article submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

31

Caiazzo, A. (2005). Analysis of Lattice Boltzmann Initialization Routines. Journal
of Statistical Physics, 121(1-2):37–48.
Chen, Q., Jiang, L., Li, C., Hu, D., Bu, J.-w., Cai, D., and Du, J.-l. (2012).
Haemodynamics-driven developmental pruning of brain vasculature in zebrafish.
PLoS Biology, 10(8):e1001374.
Chen, S. and Doolen, G. D. (1998). Lattice Boltzmann method for fluid flows.
Annual Review of Fluid Mechanics, 30(1):329–364.
Chien, S. (1970). Shear Dependence of Effective Cell Volume as a Determinant of
Blood Viscosity. Science, 168(3934):977–979.
Feke, G. T., Tagawa, H., Deupree, D. M., Goger, D. G., Sebag, J., and Weiter, J. J.
(1989). Blood flow in the normal human retina. Investigative ophthalmology &
visual science, 30(1):58–65.
Formaggia, L., Quarteroni, A., and Veneziani, A. (2009). Cardiovascular Mathematics. Modeling and simulation of the circulatory system. Springer-Verlag Italia,
Milano.
Franco, C. A., Mericskay, M., Parlakian, A., Gary-Bobo, G., Gao-Li, J., Paulin, D.,
Gustafsson, E., and Li, Z. (2008). Serum response factor is required for sprouting
angiogenesis and vascular integrity. Developmental cell, 15(3):448–61.
Fruttiger, M. (2007).
10(2):77–88.

Development of the retinal vasculature.

Angiogenesis,

Ganesan, P., He, S., and Xu, H. (2010). Development of an image-based network
model of retinal vasculature. Annals of biomedical engineering, 38(4):1566–85.
Gariano, R. F. and Gardner, T. W. (2005). Retinal angiogenesis in development
and disease. Nature, 438(7070):960–6.
Geudens, I. and Gerhardt, H. (2011). Coordinating cell behaviour during blood
vessel formation. Development, 138(21):4569–83.
Groen, D., Hetherington, J., Carver, H. B., Nash, R. W., Bernabeu, M. O., and
Coveney, P. V. (2013). Analysing and modelling the performance of the HemeLB
lattice-Boltzmann simulation environment. Journal of Computational Science,
4(5):412–422.
Hardy, P., Abran, D., Li, D. Y., Fernandez, H., Varma, D. R., and Chemtob, S.
(1994). Free radicals in retinal and choroidal blood flow autoregulation in the
piglet: interaction with prostaglandins. Investigative ophthalmology & visual science, 35(2):580–91.
Heywood, J. G., Rannacher, R., and Turek, S. (1996). Artificial boundaries and
flux and pressure conditions for the incompressible Navier-Stokes equations. International Journal for Numerical Methods in Fluids, 22(5):325–352.
Jones, E. a. V., le Noble, F., and Eichmann, A. (2006). What determines blood
vessel structure? Genetic prespecification vs. hemodynamics. Physiology, 21:388–
95.
Article submitted to Royal Society

32

Bernabeu, Jones, Nielsen, et al.

Kiel, J. W. (2010). The Ocular Circulation. Morgan & Claypool Life Sciences.
Krüger, T., Varnik, F., and Raabe, D. (2009). Shear stress in lattice Boltzmann
simulations. Physical Review E, 79(4):046704.
Ladd, A. and Verberg, R. (2001). Lattice-Boltzmann simulations of particle-fluid
suspensions. Journal of Statistical Physics, 104(September):1191–1251.
Ladd, A. J. C. (1994). Numerical simulations of particulate suspensions via a
discretized Boltzmann equation. Part 1. Theoretical foundation. Journal of Fluid
Mechanics, 271:285.
Lätt, J. and Chopard, B. (2008). Straight velocity boundaries in the lattice Boltzmann method. Physical Review E, 77(5):056703.
Lawson, N. D. and Weinstein, B. M. (2002). Arteries and veins: making a difference
with zebrafish. Nature reviews. Genetics, 3(9):674–82.
Liu, D., Wood, N. B., Witt, N., Hughes, A. D., Thom, S. A., and Xu, X. Y. (2009).
Computational Analysis of Oxygen Transport in the Retinal Arterial Network.
Current Eye Research, 34(11):945–956.
Lobov, I. B., Rao, S., Carroll, T. J., Vallance, J. E., Ito, M., Ondr, J. K., Kurup, S.,
Glass, D. a., Patel, M. S., Shu, W., Morrisey, E. E., McMahon, A. P., Karsenty,
G., and Lang, R. a. (2005). WNT7b mediates macrophage-induced programmed
cell death in patterning of the vasculature. Nature, 437(7057):417–21.
Mazzeo, M. D. and Coveney, P. V. (2008). HemeLB: A high performance parallel
lattice-Boltzmann code for large scale fluid flow in complex geometries. Computer
Physics Communications, 178(12):894–914.
Mei, R., Luo, L.-S., Lallemand, P., and D’Humières, D. (2006). Consistent initial
conditions for lattice Boltzmann simulations. Computers & Fluids, 35(8-9):855–
862.
Nagaoka, T. and Yoshida, A. (2006). Noninvasive evaluation of wall shear stress on
retinal microcirculation in humans. Investigative ophthalmology & visual science,
47(3):1113–9.
Nash, R. W., Carver, H. B., Bernabeu, M. O., Hetherington, J., Groen, D., Krüger,
T., and Coveney, P. V. (2014). Choice of boundary condition and collision operator for lattice-Boltzmann simulation of moderate Reynolds number flow in
complex domains. Physical Review E, 89:023303.
Ninomiya, H. and Inomata, T. (2006). Microvasculature of the mouse eye: Scanning
electron microscopy of vascular corrosion casts. Journal of Experimental Animal
Science, 43(3):149–159.
Paques, M. (2003). Structural and Hemodynamic Analysis of the Mouse Retinal
Microcirculation. Investigative Ophthalmology & Visual Science, 44(11):4960–
4967.
Article submitted to Royal Society

Computational analysis of haemodynamics during angiogenesis

33

Popel, A. S. and Johnson, P. C. (2005). Microcirculation and Hemorheology. Annual
review of fluid mechanics, 37:43–69.
Potente, M., Gerhardt, H., and Carmeliet, P. (2011). Basic and therapeutic aspects
of angiogenesis. Cell, 146(6):873–87.
Pries, A. R., Secomb, T. W., and Gaehtgens, P. (1996). Biophysical aspects of
blood flow in the microvasculature. Cardiovascular research, 32(4):654–67.
Schroeder, W. J., Martin, K., and Lorensen, W. E. (2003). The Visualization
Toolkit: An Object-Oriented Approach to 3D Graphics, Third Edition.
Siekmann, A. F. and Lawson, N. D. (2007). Notch signalling limits angiogenic cell
behaviour in developing zebrafish arteries. Nature, 445(7129):781–4.
Siggers, J. H. and Waters, S. L. (2005). Steady flows in pipes with finite curvature.
Physics of Fluids, 17(7):077102.
Stahl, B., Chopard, B., and Latt, J. (2010). Measurements of wall shear stress
with the lattice Boltzmann method and staircase approximation of boundaries.
Computers & Fluids, 39(9):1625–1633.
Uemura, A., Kusuhara, S., Katsuta, H., and Nishikawa, S.-I. (2006). Angiogenesis
in the mouse retina: a model system for experimental manipulation. Experimental
cell research, 312(5):676–83.
Vogel, J., Kiessling, I., Heinicke, K., Stallmach, T., Ossent, P., Vogel, O., Aulmann, M., Frietsch, T., Schmid-Schönbein, H., Kuschinsky, W., and Gassmann,
M. (2003). Transgenic mice overexpressing erythropoietin adapt to excessive
erythrocytosis by regulating blood viscosity. Blood, 102(6):2278–84.
Wang, Y., Bower, B. a., Izatt, J. a., Tan, O., and Huang, D. (2007). In vivo total
retinal blood flow measurement by Fourier domain Doppler optical coherence
tomography. Journal of biomedical optics, 12(4):041215.
Wang, Y., Fawzi, A., Tan, O., Gil-Flamer, J., and Huang, D. (2009). Retinal blood
flow detection in diabetic patients by Doppler Fourier domain optical coherence
tomography. Optics Express, 17(5):4061.
Watson, M. G., McDougall, S. R., Chaplain, M. a. J., Devlin, a. H., and Mitchell,
C. a. (2012). Dynamics of angiogenesis during murine retinal development: a
coupled in vivo and in silico study. Journal of the Royal Society, Interface / the
Royal Society, 9(74):2351–64.
White, B., Pierce, M., Nassif, N., Cense, B., Park, B., Tearney, G., Bouma, B.,
Chen, T., and de Boer, J. (2003). In vivo dynamic human retinal blood flow
imaging using ultra-high-speed spectral domain optical coherence tomography.
Optics express, 11(25):3490–7.
Windberger, U., Bartholovitsch, A., Plasenzotti, R., Korak, K., and Heinze, G.
(2003). Whole blood viscosity, plasma viscosity and erythrocyte aggregation in
nine mammalian species: reference values and comparison of data. Exp Physiol,
88(3):431–440.
Article submitted to Royal Society

34

Bernabeu, Jones, Nielsen, et al.

Windberger, U., Grohmann, K., Goll, a., Plasenzotti, R., and Losert, U. (2005).
Fetal and juvenile animal hemorheology. Clinical hemorheology and microcirculation, 32(3):191–7.
Wright, W. S. and Harris, N. R. (2008). Ozagrel attenuates early streptozotocininduced constriction of arterioles in the mouse retina. Experimental eye research,
86(3):528–36.
Wright, W. S., Messina, J. E., and Harris, N. R. (2009). Attenuation of diabetesinduced retinal vasoconstriction by a thromboxane receptor antagonist. Experimental eye research, 88(1):106–12.
Wright, W. S., Yadav, A. S., McElhatten, R. M., and Harris, N. R. (2012). Retinal
blood flow abnormalities following six months of hyperglycemia in the Ins2(Akita)
mouse. Experimental eye research, 98:9–15.
Xiong, W. and Zhang, J. (2010). Shear stress variation induced by red blood cell
motion in microvessel. Annals of biomedical engineering, 38(8):2649–59.
Yazdanfar, S. (2003). In Vivo Imaging of Human Retinal Flow Dynamics by Color
Doppler Optical Coherence Tomography. Archives of Ophthalmology, 121(2):235.
Zhi, Z., Cepurna, W., Johnson, E., Shen, T., Morrison, J., and Wang, R. K. (2011).
Volumetric and quantitative imaging of retinal blood flow in rats with optical
microangiography. Biomedical optics express, 2(3):579–91.
Zhi, Z., Yin, X., Dziennis, S., Wietecha, T., Hudkins, K. L., Alpers, C. E., and Wang,
R. K. (2012). Optical microangiography of retina and choroid and measurement
of total retinal blood flow in mice. Biomedical optics express, 3(11):2976–86.
Ziegler, D. P. (1993). Boundary conditions for lattice Boltzmann simulations. Journal of Statistical Physics, 71(5-6):1171–1177.

Article submitted to Royal Society

