The distribution of atherosclerotic lesions within the rabbit vasculature, particularly within the descending thoracic aorta, has been mapped in numerous studies. The patchy nature of such lesions has been attributed to local variation in the pattern of blood flow. However, there have been few attempts to model and characterize the flow. In this study, a high-order continuous Galerkin finite-element method was used to simulate blood flow within a realistic representation of the rabbit aortic arch and descending thoracic aorta. The geometry, which was obtained from computed tomography of a resin corrosion cast, included all vessels originating from the aortic arch (followed to at least their second generation) and five pairs of intercostal arteries originating from the proximal descending thoracic aorta. The simulations showed that small geometrical undulations associated with the ductus arteriosus scar cause significant deviations in wall shear stress (WSS). This finding highlights the importance of geometrical accuracy when analysing WSS or related metrics. It was also observed that two Dean-type vortices form in the aortic arch and propagate down the descending thoracic aorta (along with an associated skewed axial velocity profile). This leads to the occurrence of axial streaks in WSS, similar in nature to the axial streaks of lipid deposition found in the descending aorta of cholesterol-fed rabbits. Finally, it was observed that WSS patterns within the vicinity of intercostal branch ostia depend not only on local flow features caused by the branches themselves, but also on larger-scale flow features within the descending aorta, which vary between branches at different locations. This result implies that disease and WSS patterns in the vicinity of intercostal ostia are best compared on a branch-by-branch basis.
The prevalence of atherosclerosis varies within the mammalian vasculature, particularly in regions of arterial branching and curvature. Since these regions are associated with complex blood flow patterns, it has been postulated that blood flow may play an important role in regulating atherogenesis. Fry  inferred from studies of hypercholesterolaemic animals that lesions develop in regions of high wall shear stress (WSS), hypothesizing that high WSS damaged the arterial endothelium, and hence increased lipoprotein influx into the intima. However, Caro et al. [2,3], who studied human post-mortem material, suggested that lesions develop in regions of low wall shear rate (WSR), owing to a decreased efflux of material made or modified within the arterial wall.
The low WSR theory of Caro et al. is widely accepted today (although the original mechanism is contradicted by evidence that lesion prevalence depends on plasma low-density lipoprotein (LDL) concentrations). However, the existence of a direct correlation between disease patterns and WSR/WSS magnitude may be too simplistic. More complex relationships have been suggested, including the proposal that a combination of low and temporally oscillating WSR/WSS predisposes to disease [4–6]. A better understanding of arterial haemodynamics is required before definitive correlations between flow and disease can be established. In recent years, understanding has been enhanced by numerous computational studies [7–12]. Such studies circumvent technical limitations associated with measuring flow in vivo and facilitate the analysis of derived metrics such as WSR/WSS. However, despite numerous experimental studies of disease location in the rabbit vasculature, there have been few computational studies of blood flow in this species and the corresponding flow therefore remains largely uncharacterized.
In the present study, a high-order continuous Galerkin finite-element method was used to model blood flow within a realistic representation of the rabbit aortic arch and descending thoracic aorta. The geometry was obtained from high-resolution computed tomography (CT) of a vascular cast. It included all vessels originating from the aortic arch (followed to at least their second generation) and five pairs of intercostal arteries originating from the proximal descending thoracic aorta. Steady-state flow of a Newtonian fluid was modelled within the geometry, which was assumed to have rigid and stationary walls. These assumptions, which considerably simplified calculations, clearly represent a significant approximation of in vivo conditions. However, they permitted understanding of the underlying steady-state fluid dynamics, and did not prevent important new information being gained owing to the unprecedented detail and extent of the geometrical definition.
2. Method: obtaining a geometrical definition of the rabbit aorta
The following section details the process of obtaining a realistic geometrical definition and volume mesh of the rabbit aortic arch and descending thoracic aorta. The process consists of four stages: casting the rabbit vasculature, CT imaging of the resulting cast, computational reconstruction of the CT data and computational generation of a curved element volume mesh from the surface definition. Similar approaches have been employed in various previous studies [7,11]. However, to the authors' knowledge, no previous studies have obtained geometrical definitions with the combined detail and extent of the one presented here.
2.2. Casting procedure
A male New Zealand white rabbit aged 18 months and weighing 3.25 kg received heparin (2000 USP units, IV) and was then euthanized with an overdose of sodium pentobarbitone (Euthatal, 160 mg kg−1, IV). A cannula was inserted and secured into the cardiac left ventricle. It was used to flush the arterial system with saline and then to infuse Batsons No. 17 casting material (Polysciences, Inc.) at a pressure of 95–100 mm Hg. This pressure, which is physiologically realistic for the rabbit , was maintained until the casting material hardened and could no longer be injected (at which point approx. 70 ml had been used). During the initial exothermic stage of the curing process, cold water was trickled through the chest and abdominal cavities to dissipate heat. The casting material was then left overnight to cure fully.
The following day, the rabbit carcass was submerged in a potassium hydroxide solution (approx. 30% w/v). After two weeks, the solution was renewed and left for a further 24 h. The second potassium hydroxide solution was then replaced with a detergent solution and left for another 24 h. Finally, the cast was removed from the detergent and gently washed with warm water. Images of the cast are shown in figure 1. As an indication of scale, the aortic root has a diameter of approximately 6 mm.
2.3. Imaging procedure
The cast was imaged using a Metris X-Tek HMX-ST Micro-CT Scanner at the Natural History Museum, London, UK. A 2000 × 2000 × 2000 voxel cube of DICOM data (approx. 9.2 GB) was obtained, with a voxel size of 55.2 µm. Over 100 voxels spanned the 6 mm diameter of the aortic root.
2.4. Obtaining a surface definition
An intensity isosurface was extracted from the DICOM data using Amira v. 5.2.0 (Visage Imaging, Inc.) The surface was cropped using the Vascular Modelling Tool Kit (VMTK) to remove unwanted arteries. Flow extensions terminating with a circular cross section were then extruded from each outflow vessel using VMTK, to facilitate the application of boundary conditions. Finally, the surface was smoothed using VMTK. The degree of smoothing was highly limited in order to avoid loss of small-scale geometrical features. Images of the final surface definition, including flow extensions, are shown in figure 2. A VTP file of the final surface definition is available online as electronic supplementary material.
2.5. Obtaining a volume mesh
Use of a high-order continuous Galerkin finite-element method requires that the computational domain is subdivided into a mesh of elements, within which the solution is represented by a high-order polynomial expansion. Since each element holds multiple degrees of freedom, the volume mesh can be coarse compared with a typical low-order finite element or finite volume mesh. In this study, a volume mesh of 79 847 tetrahedral elements and 25 498 prismatic elements was constructed using Gambit v. 2.4.6 and Tgrid v. 4.0.24 (ANSYS, Inc.). The size of the elements was based on the local curvature of the surface definition. The prismatic elements were placed adjacent to the luminal surface of the arterial wall in order to improve resolution of any momentum boundary layers that formed. External faces of prismatic boundary elements were curved using SPHERIGON patches , following the approach of Grinberg et al. , so that they accurately approximated the surface definition despite their large size. Sections of the resulting volume mesh are shown in figure 3.
3. Method: simulating blood flow within the rabbit aorta
3.1. Governing equations
In the present study, blood was treated as an incompressible Newtonian fluid. Also, time-dependent dynamics were ignored, and only steady-state simulation were performed. Specifically, blood flow was modelled using the incompressible steady-state Navier–Stokes equations for a Newtonian fluid, which can be written as 3.1and 3.2where μ is the viscosity of rabbit blood, ρ is the density of rabbit blood, u is the three-dimensional velocity vector field within the rabbit aorta and p is the pressure field within the rabbit aorta. Such an approach has been adopted in numerous previous studies [10,12].
There are various reasons why steady-state simulations were performed, as opposed to fully time-dependent studies. Primarily, the results of steady-state simulations can be analysed with relative ease; giving insight that will facilitate interpretation of future (more complex) unsteady results, and focus the remit of future (more costly) unsteady studies. Additionally, it was felt that insufficient experimental data existed to reliably parameterize time-dependent conditions at each of the 22 open boundaries (1 inlet and 21 outlets). In reality, however, blood flow is pulsatile. Consequently, there are disadvantages associated with steady-state simulations. In particular, they preclude calculation and analysis of potentially important metrics such as oscillatory shear index (OSI) [4,16]. Also, they cannot capture truly unsteady flow phenomena (but only quasi-steady dynamics). The potential effect of pulsatility on the main findings of the present study will be discussed towards the end of the paper in §5.
3.2. Boundary conditions
A spatially uniform inflow velocity profile was prescribed at the aortic root. Applying such a profile is an approximation, since flow from the left ventricle is likely to be skewed and possess an in-plane component. Parabolic (fully developed) outflow profiles were prescribed perpendicular to all other boundaries, apart from the main outflow of the descending thoracic aorta, where a zero velocity gradient condition was prescribed perpendicular to the outflow plane. At the arterial wall, which was assumed to be rigid and stationary, a no-slip boundary condition was applied (i.e. all components of the velocity were set to zero).
3.3. Parameter values
3.3.1. Reynolds numbers
For any viscous flow, one can define a Reynolds number (Re), which quantifies the ratio of inertial to viscous forces acting within the fluid. In this study, Re is defined as 3.3where D0 is the hydraulic diameter of the aortic root inflow plane (the diameter of a circle with the same area as the aortic root inflow plane), and V0 is the spatially constant magnitude of the inflow velocity at the aortic root.
The average flow rate into the rabbit ascending thoracic aorta has been measured by Avolio et al.  to be 3.3 ml s−1. These measurements were made using an electromagnetic probe with an internal diameter of 5 mm that caused ‘minimal constriction’ of the ascending aorta. Assuming that the arterial wall thickness was 7 per cent of the vessel diameter (this is constant across many species and vessels ), then the inflow area was approximately 15 mm2 and the measured inflow rate corresponds to an inflow velocity magnitude of approximately 0.22 m s−1. The peak inflow velocity magnitude measured by Avolio et al.  was 0.927 m s−1. In the present study, three values of inflow velocity magnitude V0 were considered. Specifically, V0 = 0.0332 m s−1 (an arbitrary low value), V0 = 0.199 m s−1 (based on the mean flow rate measurements of Avolio et al. ) and V0 = 0.862 m s−1 (based on the peak flow rate measurements of Avolio et al. ). Assuming that ρ = 1.044 × 103 kg m−3 , μ = 4.043 × 10−3 kg m−1 s−1  and D0 = 5.838 × 10−3 m (measured from the geometry), these values of V0 lead to round values Re = 50, Re = 300 and Re = 1300, respectively.
3.3.2. Flow splits
Flow splits into the innominate artery and the left subclavian artery were based on the measurements of Barakat et al. . Specifically, 14.7 per cent of the inflow at the aortic root went into the innominate artery and 7.1 per cent into the left subclavian artery. Flow splits for all subsequent branches of these arteries were determined using Murray's law , which states that the flow in a vessel is proportional to the diameter of the vessel cubed. It was assumed that on average, each intercostal branch received 0.2 per cent of the flow entering the descending thoracic aorta. This value sits within the range considered by Kazakidi et al.  but is otherwise arbitrary. The 10 intercostal branches, therefore, received between them a total of 2 per cent of the flow entering the descending thoracic aorta (which is equivalent to 1.56% of the flow entering at the aortic root). The exact split of this aggregate flow between the individual intercostal branches was once again determined using Murray's law.
3.4. Numerical method
The incompressible steady-state Navier–Stokes equations were solved within the arterial geometry using a high-order continuous Galerkin method. This approach combines the superior accuracy of spectral methods with the geometrical flexibility of low-order finite element or finite volume schemes. For further details of the method, see papers by Sherwin & Ainsworth  and Sherwin & Karniadakis , the textbook by Karniadakis & Sherwin  and references therein.
Convergence of the solution was verified by performing simulations at Re = 300 using third-, fifth- and seventh-order polynomial expansions to represent the solution within each element. Since the computational domain contains 79 847 tetrahedral elements and 25 498 prismatic elements, use of these expansions results in a total of approximately 3.07 × 106, 1.37 × 107 and 3.71 × 107 degrees of freedom, respectively (accounting for vertex, edge and face multiplicity, and counting each of the velocity components and the pressure as separate degrees of freedom). When switching from third-order to fifth-order solution polynomials, the total wall force changed by 0.019 per cent, and when switching from fifth-order to seventh-order solution polynomials, the total wall force changed by 0.0040 per cent. The solution was, therefore, deemed to have converged when using fifth-order polynomials within each element, and thus all presented results were obtained using this order. For comparison with a more standard low-order calculation, this corresponds (in terms of degrees of freedom) to using a mesh of approximately 2.05 × 107 linear tetrahedral elements. However, such as scheme will only be second-order accurate in space, whereas the scheme employed here is nominally sixth-order accurate.
4. Results and analysis
4.1. Flow in the aortic arch
Figure 4a,c shows streamlines in the aortic arch when Re = 300. Significant secondary flow was observed. The streamlines are qualitatively comparable with those obtained by Barakat et al.  using microcinematographic methods in excised rabbit aortas. Figure 4b,d shows coherent vortical structures in the aortic arch when Re = 300. The structures were defined by isosurfaces of the widely used metric λ2 , which accurately identifies vortices at low Reynolds numbers (unlike minimum-pressure criteria). Two Dean-type vortices  of unequal magnitude formed. This is consistent with literature on flow in curved pipes [29,30], and with studies of flow in the aortic arch of mice  and human subjects . However, to the knowledge of the authors, formation of such vortices has not previously been demonstrated numerically within the aortic arch of rabbits. Results for Re = 1300 (not shown) were qualitatively similar. However, at Re = 50 (results not shown), there was minimal secondary flow within the arch and no significant Dean vortices developed. In all cases, a Dean number De defined as 4.1can be associated with the vortical structures in the aortic arch, where R is the radius of curvature of the aortic arch. Measurements made using VMTK suggest R ≈ 2D0, and hence De ≈ 0.5Re. Therefore, the Reynolds numbers of Re = 50, Re = 300 and Re = 1300 correspond to Dean numbers of De = 25, De = 150 and De = 650, respectively.
Figures 5 and 6 show maps of WSS within the aortic arch for all three values of Re considered. The patterns are complex and become increasingly heterogeneous as Re is increased. Two important observations can be made. The first is that, owing to the complexity of the patterns, it would be incorrect to characterize the inner curvature of the arch as a ‘low shear region’ and the outer curvature as a ‘high shear region’ (as is common when interpreting experimental data). The second is that small surface undulations associated with the ductus arteriosus scar, while having a negligible effect on flow, induce significant deviations in WSS, comparable to those caused by vascular scale flow features. The deviations are marked with arrows in figures 5 and 6. This observation highlights the critical importance of geometrical accuracy when calculating and analysing WSS or related metrics. Interestingly, Duff  suggests that the ductus arteriosus scar is a site of atherosclerotic disease in humans.
4.2. Flow in the proximal descending thoracic aorta
Figure 7 shows plots of axial velocity and vorticity magnitude within the descending thoracic aorta when Re = 300. The axial velocity profile is skewed to one side, and the two Dean-type vortices emerging from the aortic arch continue to propagate down the descending thoracic aorta. Figure 8 shows streamlines in the descending thoracic aorta for all three values of Re considered. The results were obtained using the particle-tracking approach of Copella et al. , which was devised specifically for post-processing the results of unstructured high-order finite-element simulations. When Re = 50, there is almost no secondary flow within the descending thoracic aorta, since significant Dean-type vortices do not form in the aortic arch. However, when Re = 300 and Re = 1300, vortices are clearly present. Also, it is apparent from figure 8 that the streamline colours, which were uniformly distributed at the aortic root, become increasingly mixed as Re is increased. Hence, raising Re increases the degree of stirring to which the blood is subjected as it travels through the aortic arch and down the descending thoracic aorta.
Figure 9 shows maps of WSS in the descending thoracic aorta for all three values of Re considered. When Re = 300 and Re = 1300, the WSS patterns exhibit streaks owing to the presence of a skewed axial velocity profile. Figure 10 shows en face flattened three-colour maps of WSS in the descending thoracic aorta. In all cases, the aorta was unwrapped via a ‘computational incision’ along the ventral aortic wall, and the colour maps are presented with the luminal surface facing upwards. The patterns of WSS in figure 10 can be compared with the patterns of lipid deposition observed in the descending aorta of cholesterol-fed rabbits (figure 11). While detailed quantitative comparisons between figures 10 and 11 are unwarranted at this stage (owing to the small sample sizes), it is clear that distinctive axial streaks of both WSS and lipid deposition are present, indicating that a relationship between WSS and disease may exist. Systematic comparison of several computational and experimental results is necessary to ascertain in more detail the nature of any such relationship.
4.3. Flow into the intercostal arteries
Figure 12 shows plots of WSS in the vicinity of various intercostal branch ostia for all three values of Re considered. The patterns reflect a superposition of both ‘local’ and ‘global’ contributions. The local contributions are caused by the deflection of aortic flow into the branches themselves. Such contributions were studied in detail by Kazakidi et al.  within idealized geometries, and are expected to be approximately independent of branch location. The global contributions result from bulk flow within the descending thoracic aorta. When Re = 300 and Re = 1300, these global contributions exhibit streaks of high and low WSS, as discussed above.
When Re = 50, the WSS patterns within the vicinity of each branch are approximately independent of the branch location. This is because local contributions to the WSS are approximately independent of branch location and there is minimal spatial variation in the global contribution. However, when Re = 300 and Re = 1300, the WSS patterns within the vicinity of each branch become highly dependent on branch location, since the global contribution now varies spatially. Note that as well as varying between branch pairs, the global contributions to WSS also vary on the scale of the branch pairs themselves. This can be seen by comparing figure 12d and figure 12f which show the second and fourth pair of intercostal branches, respectively, when Re = 300. In figure 12d, one of the branches sits in a streak of high WSS, whereas the other sits on the edge of a low WSS streak. However, in figure 12f, a streak of high WSS runs across the entrance to both branches.
These findings are significant for three reasons. Firstly, they challenge the accuracy of results obtained using models that do not include the spatially varying, large-scale aortic flow field . Secondly, they may explain why disease patterns in the vicinity of intercostal branch entrances vary between different branch locations. Thirdly, they call into question the process of combining disease patterns from intercostal branches in different locations to obtain average maps . It may be more appropriate to analyse disease patterns, and compare them with predicted WSS patterns, on a branch-by-branch basis.
The main findings of the present study can be summarized as follows:
Small surface undulations associated with physiologically realistic features such as the ductus arteriosus scar (which was only resolved owing to the unprecedented detail of the geometrical definition) induce deviations in WSS comparable to those caused by vascular scale flow features. Hence, geometrical accuracy is critically important when calculating and analysing WSS or related metrics.
Dean-type vortices form in the aortic arch of rabbits. This phenomenon is consistent with studies of flow in curved pipes [29,30], and with studies of flow in the aortic arch of mice  and human subjects . However, to the knowledge of the authors, formation of such vortices has not previously been demonstrated numerically within the aortic arch of rabbits. It is observed that these Dean vortices (along with an associated skewed axial velocity profile) propagate down the descending thoracic aorta, consequently leading to the formation of axial streaks in WSS, which are similar in nature to the streaks of lipid deposition observed in the descending aorta of cholesterol-fed rabbits.
WSS patterns within the vicinity of intercostal branch ostia depend not only on local flow features caused by the branches themselves but also on larger-scale flow features within the descending aorta, which vary between branches at different locations. This result (which was only obtained owing to the combined detail and extent of the geometrical definition) implies that patterns of WSS and disease near intercostal ostia are best compared on a branch-by-branch basis.
5. The effect of pulsatility
The importance of pulsatile dynamics within the rabbit aorta can be assessed by calculating a Womersley number α, defined as 5.1where T is a period associated with the rabbit cardiac cycle. The magnitude of α represents the ratio of temporally varying inertial forces to viscous forces within the flow. If α < 1, it can be assumed that the flow is quasi-steady i.e. steady-state results will faithfully represent flow patterns at various points in the cardiac cycle. However, if α > 1, the flow may not be entirely quasi-steady. Barakat et al.  suggest that the rabbit heart rate is typically 250–300 beats per minute (making reference to the studies of Borkowski et al. , Wyatt et al. , and also the textbook of Weisbroth et al. ). Also, Avolio et al.  measured heart rates of 210 beats per minute in anaesthetized rabbits, and 288 beats per minute in unanaesthetized rabbits. Assuming a heart rate of 275 beats per minute (with all other parameter values as detailed above), a Womersley number of α ≈ 8 is obtained. Consequently, flow in the rabbit aorta is unlikely to be entirely quasi-steady.
Having established α ≈ 8, it is pertinent to consider how pulsatility (which was ignored in the present study) will affect the main findings detailed in §4.4. Finding (1) is a statement about the sensitivity of WSS to geometry, and can reasonably be considered independent of whether WSS varies in time (even if blood flow is not entirely quasi-steady). Consequently, pulsatility is unlikely to affect finding (1). The results of Pitt  can be used to estimate how pulsatility will likely affect finding (2). Pitt studied steady state and pulsatile flow in a double bend configuration at a Reynolds number of 250. In his steady-state study, Pitt observed the formation of Dean vortices, which evolved in a quasi-steady fashion when the flow pulsated with a Womersley number of four. However, when the flow pulsated with a Womersley number of eight, the vortices, while not entirely destroyed, were no longer quasi-steady. These results suggest that pulsatility may modify the Dean vortices described in finding (2), but it is unlikely to preclude their formation altogether. Additionally, in the double bend configuration described above, Pitt found that the steady-state WSS pattern was similar to the time-averaged WSS patterns for the pulsatile cases. This suggests that the streaks of WSS described in finding (2) will also occur (in a time-averaged sense at least) if the flow is pulsatile. Finally, it is reasonable to suppose that local and global WSS patterns will superimpose in the vicinity of intercostal ostia, irrespective of whether or not the flow is pulsatile. Therefore, since pulsatility is unlikely to homogenize WSS in the descending thoracic aorta, it can be concluded that pulsatility is unlikely to affect finding (3).
The analysis presented above suggests that the findings of the present study will not be significantly affected by pulsatility. However, future time-dependent studies are required to confirm this is the case. In addition, such time-dependent studies would allow characterization of the potentially non-quasi-steady Dean vortex dynamics in the descending aorta, as well as the distribution of potentially important time-dependent flow metrics such as OSI [4,16].
A high-order continuous Galerkin finite-element method was used to simulate blood flow within a realistic representation of the rabbit aortic arch and descending thoracic aorta. Steady-state flow of a Newtonian fluid was modelled within the geometry, which was assumed to have rigid and stationary walls. These assumptions clearly represent a significant simplification of in vivo conditions. Nevertheless, the detail and the extent of the geometrical definition allowed new and interesting results to be obtained. Firstly, small geometrical undulations associated with the ductus arteriosus scar caused significant deviations in WSS, highlighting the importance of geometrical accuracy when analysing WSS or related metrics. Secondly, two Dean-type vortices form in the aortic arch and (along with an associated skewed axial velocity profile) propagate down the descending thoracic aorta. This leads to the formation of axial streaks in WSS, which are similar in nature to the streaks of early-stage atherosclerosis observed to form in the descending aorta of cholesterol-fed rabbits. Thirdly, WSS patterns within the vicinity of intercostal ostia depend not only on local flow features (caused by the branches themselves), but also on larger-scale global flow features within the descending aorta (which vary between ostia at different locations). This finding has several implications. It challenges the accuracy of results obtained using models that do not include the spatially varying large-scale aortic flow field . It may also explain why disease patterns in the vicinity of intercostal ostia vary between branch locations. Finally, it calls into question the process of combining disease patterns from branches at different locations in order to obtain average maps ; it may be more appropriate to analyse disease patterns, and compare them with predicted WSS patterns, on a branch-by-branch basis.
Future studies should build on the work presented in this paper to develop a detailed and holistic understanding of blood flow within the rabbit aortic arch and descending thoracic aorta. Such an understanding is required before definitive correlations between blood flow patterns and atherosclerosis can be established with certainty. In particular, future studies should investigate the effect of varying flow splits between the aorta and its daughter vessels, the effect of using a more realistic inflow velocity profile at the aortic root, the effect of anatomical variation between animals and the effect of including more complex physical phenomena, such as pulsatile flow, deformable walls, non-Newtonian rheology and the transport of dissolved solutes.
All procedures complied with the UK Animals (Scientific Procedures) Act 1986 and were approved by the Imperial College London local ethical review panel.
This study was funded by the Engineering and Physical Sciences Research Council, the Biotechnology and Biological Sciences Research Council, Fundación Caja Madrid and a British Heart Foundation Research Excellence Award. The authors would like to thank Dr Richie Abel and Dr Stig Walsh for imaging the cast, Dr Donal Taylor and Dr Leopold Grinberg for assistance obtaining the high-order volume mesh, Miss Véronique Peiffer for assistance obtaining the unwrapped shear map images, Dr Stephanie Cremers for providing the disease images, Prof. Kim Parker for useful discussions, and the high-performance computing staff at Imperial College London for useful advice. Finally, the authors would like to thank the Barcelona Supercomputing Centre (HPC-EUROPA2 project number 228398), and Prof. George Karniadakis, for use of their facilities.
- Received March 2, 2011.
- Accepted April 21, 2011.
- This journal is © 2011 The Royal Society