## Abstract

High-resolution unsteady three-dimensional flow simulations in large intracranial arterial networks of a healthy subject and a patient with hydrocephalus have been performed. The large size of the computational domains requires the use of thousands of computer processors and solution of the flow equations with approximately one billion degrees of freedom. We have developed and implemented a two-level domain decomposition method, and a new type of outflow boundary condition to control flow rates at tens of terminal vessels of the arterial network. In this paper, we demonstrate the flow patterns in the normal and abnormal intracranial arterial networks using patient-specific data.

## 1. Introduction

The dynamics of blood flow in the human brain depend upon a complex network of vessels under a variety of temporal and spatial constraints. Abnormalities in the delivery of blood to the brain clearly underlie the pathophysiology of stroke, vasospasm, traumatic brain injury, vascular dementias and probably conditions such as migraine and hydrocephalus. Clinical decisions are often made on the basis of steady-state conditions (e.g. mean intracranial pressures, mean cerebral blood flow, etc.), but there is clearly a risk that ignoring the range of spatial and temporal scales present may limit understanding, and hence clinical effectiveness. Considerable interest attends the recent observation that a frequency-dependent pulsation absorber may be a part of normal physiology, dampening the effect of cardiac pulsation on microvessels (Zou *et al*. 2008). A detailed understanding of such a phenomenon categorically depends on the development of new techniques to grasp the importance of dynamic processes at a variety of different scales.

The range of length and times scales that needs to be considered spans several orders of magnitude, from flow in large arteries to subcellular processes. Interactions of blood flow in the human body occur between different scales, in which the large-scale flow features are being coupled to cellular and subcellular biology, or at similar scales in different regions of the vascular system. The human brain is a fascinating object of study in this regard since 20 per cent of cardiac output must be distributed to tissues with exquisite regulation from a pulsatile cardiac pump, with the unique boundary condition that the intracranial volume is bounded by the skull, a rigid exoskeleton. The typical anatomy includes delivery of blood by two internal carotid (ICA) and two vertebral arteries, which branch and link to form a ring—the circle of Willis (CoW)—which can potentially provide alternative supply to any area of the brain if one or more of the supply trunks is interrupted. Abnormalities in the CoW are not uncommon, affecting up to 50 per cent of the population, according to Lippert & Pabst (1985). However, the effectiveness of the normal CoW arrangement, or the consequence of an abnormality, relies on downstream responses, and thus critically on interaction between events observable at different scales. Development of tools to simulate the events *simultaneously at many scales* becomes critical to understanding, and perhaps clinical effectiveness. We have focused on hydrocephalus, where interactions between pressures and flows in large compartments (the ventricles) may have a clinical manifestation through changes in flow pattern in microvessels (Zou *et al*. 2008).

We can begin to probe differences between dynamics in the normal state and hydrocephalus by looking at normal versus abnormal flows in an unusual subject with hydrocephalus and a CoW abnormality. The goal here is to demonstrate the technical possibilities. Application of these techniques to more subtle abnormalities may clarify the true dynamic abnormality in hydrocephalus, and perhaps many other vascular conditions where changes in the topology and geometry of the vascular tree may directly impact the risk of later severe clinical events such as ischaemic stroke and haemorrhage. Examples include abnormal arteriovenous fistulae and complex connections such as arteriovenous malformations (Friedlander 2007), atherosclerotic narrowing of vessels (Amin-Hanjani *et al*. 2005) and moyamoya syndrome (Scott *et al*. 2004; Smith & Scott 2005), where a consequence of large vessel occlusion in childhood results in the development of a dense network of microvessels to replace the function of an occluded part of the vascular tree. Therapeutic goals in these cases involve occluding abnormal channels, stenting partially occluded channels and providing alternative routes for blood to revascularize the tissue, respectively. Planning for pharmacological, surgical or endovascular treatment of these abnormalities of blood vessels could all benefit from insights derived from realistic models of the circulatory dynamics in the diseased arterial network. Preliminary work on the use of quantitative magnetic resonance angiography (MRA) to predict stroke appears promising (Amin-Hanjani *et al*. 2005).

Full three-dimensional simulation of the intracranial flow and its interaction with the brain tissues at *all* length scales is not feasible, and hence a new hierarchical modelling approach for the human arterial tree is required; such modelling is characterized by three distinct spatial length scales given below.

*The macrovascular network*(MaN) consisting of large arteries, down to diameters of 0.5 mm, which are patient specific and can be reconstructed from clinical imaging. An example of a complex three-dimensional intracranial arterial network is presented in figure 1.*The mesovascular network*(MeN) consisting of small arteries and arterioles, from 10 to 500 μm, which follow a tree-like structure governed by specific fractal laws (Zamir 1999). For example, the human brain contains approximately 10 million small arteries and arterioles.*The microvascular network*(MiN) consisting of the capillary bed, which follows a net-like structure. Its topological statistics have recently been quantified for the human brain in Cassot*et al*. (2006). The typical number of capillary segments in the brain is more than a billion.

The simulations of MaN scale are required to study three-dimensional flow dynamics, while MeN and MiN scale problems serve as a closure for the MaN simulations by providing proper boundary conditions.

The challenge of simulating the blood flow and blood flow interactions, from the scale of the entire human arterial tree down to the scale of capillary bed, lies in the high demand for supercomputing. It requires development of new numerical algorithms, parallel solvers and mathematical models for all three length scales, and also methods for validation of numerical results. In this work, we concentrate on blood flow simulations at the MaN scale, and present results of two simulations: flow in complete CoW of a healthy subject and in incomplete CoW of a patient with multiple brain disorders. We use patient-specific geometrical models of the arterial networks, measured *in vivo* flow rates for the inlet boundary conditions and the resistance–capacitance (RC)-type outflow boundary condition as a closure problem for the three-dimensional simulation. We note that the major impediment to large-scale three-dimensional flow simulation is solving an extremely large problem (in terms of the number of degrees of freedom) within a reasonable amount of time. High numerical resolution is important in order to eliminate the numerical discretization errors and evaluate precisely the sensitivity of the solution to modelling parameters. Accurate and fast reconstruction of arterial geometry, efficient management of terabytes of output data and visualization of numerical solution are also non-trivial tasks.

This paper is organized as follows. In §2, we overview the existing methods for arterial flow simulations; in §3, we present the results of two brain blood flow simulations, the first corresponding to the healthy subject with a complete CoW and the second to a patient who has incomplete CoW and also suffers from hydrocephalus; and in §4, we provide a summary and discuss goals for future research. The results presented in this paper have been obtained with our own research code NEKTAR (http://www.cfm.brown.edu/crunch/nektar.html).

## 2. Methods

Blood flow simulations at similar scales in different regions of the vascular system require reconstruction of very large and complex arterial networks, counting hundreds or even thousands of arteries. To simulate unsteady three-dimensional flow in such networks, considerable computational resources are required, lack of which in the past has led to developing reduced models of the blood circulation. Traditionally, flow interaction in large arterial trees was performed using computationally inexpensive one-dimensional models (Avolio 1980; Peiro *et al*. 2003; Sherwin *et al*. 2003*a*,*b*) or zero-dimensional (lumped) models (Formaggia *et al*. 2004). Although numerically efficient, the one-dimensional simulations can be applied successfully where the flow does not exhibit significant three-dimensional effects, e.g. in straight and narrow vessels. For example, the one-dimensional model cannot be applied for flow simulations in an aneurysm. The arterial wall compliance, usually incorporated in the one-dimensional models, requires knowledge of the mechanical properties of the arterial walls, such as Young's modulus and wall thickness. In patient-specific simulations, such properties are often unknown. Xiu & Sherwin (2007) used stochastic one-dimensional modelling to compensate for the uncertainty. Reduction of the three-dimensional flow equations to one-dimensional ones requires certain assumptions on the velocity profile to model energy losses due to friction. In steady pipe flow, such a profile is parabolic, unlike unsteady flow in complex geometry. Usually, quasi-stationary flow (parabolic velocity profile) or uniform flow (flat profile) is assumed.

The three-dimensional simulations are essential in capturing secondary flows, turbulence, recirculation and low-residence time zones, wall shear stress distribution, oscillatory stress index, etc. Owing to the geometric complexity of the arterial system, finite-element (FE) discretization of computational domain is preferred in solving the governing flow equations. High-order spectral/*hp* element methods, employed in this study, use similar to the FE discretization of the domain; however, inside each element, the solution is approximated with high-order polynomial expansion, which provides high spatial resolution. Use of high-order spectral/*hp* discretization is particularly valuable in simulations of transitional and turbulent flows (Fischer *et al*. 2007; Varghese *et al*. 2007). A key numerical issue in three-dimensional large-scale simulations is the use of fast parallel solvers and preconditioners scalable to thousands of processors (Sherwin & Casarin 2001; Tufo & Fischer 2001; Bergen *et al*. 2005; Grinberg *et al*. in press).

In the past, the three-dimensional flow solvers were typically employed for simulations in relatively small domains counting a few vessels; such simulations required only a few hundred processors to obtain a fast solution of the governing equation at each time step. However, several attempts to simulate blood flow dynamics in a tree of tens of arteries and bifurcations have been reported in the literature (Dong *et al*. 2006*b*; Vignon-Clementel *et al*. 2006; Grinberg *et al*. 2009).

Simulation of the flow dynamics in very large three-dimensional arterial trees requires use of tens or hundreds of thousands of computer processors. A simple estimate of the enormous resolution required can be obtained as follows. Assuming that we use tetrahedral elements to discretize the typical blood volume (5 l), with a 0.5 mm edge for each tetrahedron (and corresponding volume of 0.0147 mm^{3}), we will require more than 339 million finite elements. High-order spectral discretization requires construction of (*P*+3)(*P*+2)^{2} point grid inside each tetrahedral element using Gauss–Lobatto points in one direction and Gauss–Radau points in the other two directions (Karniadakis & Sherwin 1999). For relatively low-order polynomial expansion, *P*=4, the total number of grid points (number of unknowns per one variable) will be 85.5 billion—an enormous number indeed even for computations on the most advanced supercomputers. To significantly reduce the problem size, further simplification can be made, for example modelling the venous system as a large blood reservoir and treating the remaining arterial network using the three-level MaN–MeN–MiN approach. But even with the suggested simplifications, the size of the three-dimensional domain remains such that tens of thousands of computer processors must be employed to perform the large-scale high-resolution flow simulation in *reasonable* time. Solutions of the MeN and MiN problems require stochastic modelling due to high level of uncertainty in modelling parameters, such as arterial elasticity, effective viscosity, geometry, etc. Stochastic simulations of a flow in millions of arterioles will also demand significant computational resources.

A straightforward approach for the simulation of a three-dimensional flow on a MaN scale is computationally prohibitive even on petaflop computers, due to the extremely large size of the tightly coupled problem and high cost of *blocking* communication over thousands of processes. In order to exploit the available computational resources efficiently, new numerical algorithms must be developed of the type we present in §2*a*.

Simulation of flow in arterial networks with multiple outlets requires a careful approach for imposing outflow boundary conditions. The outflow boundary conditions typically include the Neumann boundary condition for the velocity fields and the Dirichlet boundary condition for the pressure. The pressure boundary condition affects the flow distribution in the network bifurcations. The importance of outflow boundary conditions has been long recognized by the community and the large literature reflects the same (Olufsen 1999; Sherwin *et al*. 2003*b*; Smith *et al*. 2004; Formaggia *et al*. 2006; Vignon-Clementel *et al*. 2006; Spilker *et al*. 2007). There are several methods to impose the pressure boundary conditions, and the four main approaches are: (i) constant pressure boundary condition, (ii) resistance boundary condition (Sherwin *et al*. 2003*b*; Formaggia *et al*. 2006; Vignon-Clementel *et al*. 2006), (iii) windkessel model boundary condition (Gibbons & Shadwick 1991; Nichols *et al*. 1998; Formaggia *et al*. 2006), and (iv) impedance boundary condition (Vignon-Clementel *et al*. 2006; Spilker *et al*. 2007; Steele *et al*. 2007). The resistance, windkessel and impedance boundary conditions are derived from modelling the peripheral resistance and compliance, and require knowledge of the peripheral arterial network characteristics. To perform patient-specific blood flow simulations and incorporate clinically measured flow rates, Grinberg & Karniadakis (2008*a*) suggested a new numerical method, which we overview in §2*b*.

To provide patient-specific boundary conditions, the blood flow rates in the arteries of interest were measured using two methods: (i) Doppler ultrasound and (ii) phase-contrast magnetic resonance imaging (PC-MRI). The first method is based on measuring the vessel diameter and the velocity at one or several points and relating the velocity to the flow rate. The second method is more complex, but it provides high spatial resolution and can be applied to measure flow velocity in arteries located behind thick bones, which are typically invisible to Doppler ultrasound. Application of PC-MRI is described in §2*c*.

### (a) Two-level domain decomposition

The development of the two-level domain decomposition (2DD) was motivated by advances in arterial flow simulations on distributed computers. The first simulation of blood flow in multiple arterial segments coupled with real-time visualization was performed during the supercomputing exhibition in Seattle in 2005 by Dong *et al*. (2006*a*); this work considered a loose coupling of the one- and three-dimensional models. The flow rates predicted by one-dimensional simulation were imposed as inlet boundary conditions for *discontinuous* three-dimensional domains. The solution in the three-dimensional domains was computed by hundreds of processors on four supercomputers, geographically distributed across the USA and also on one supercomputer in the UK. The communication was handled using grid technology by MPICH-g2 library (recently replaced by the more advanced MPIg library). The requirement for distributed computing was due to an insufficient number of computer processors offered by each of the supercomputing centres. The new generation of massively parallel supercomputers allows performing the simulations within a single machine. From the perspective of parallel computing, this simulation was a significant achievement (Boghosian *et al*. 2006).

In our simulations, we consider very large *continuous* three-dimensional computational domains counting tens of arteries and bifurcations. In order to perform such large-scale simulations, we have developed a new domain decomposition method (Grinberg & Karniadakis 2008*b*). The method is based on two-level partitioning: the coarse and the refined partitioning. On the coarse level, the domain is subdivided into a series of large subdomains (patches), then in each patch a standard domain partitioner is applied. Multi-level partitioning of the computational domain requires multi-level parallelism in order to maintain high parallel efficiency in simulations on thousands of processors of a single computer and, particularly, on distributed computers. Typically, we employ four to five layers of processor subgroups forming a multi-level communicating interface; the blocking communication between processors is performed within non-overlapping groups only and non-blocking communication between the subgroups is significantly minimized. Reducing the number of processors involved in blocking communications enhances the effectiveness of the parallel solver.

The entire simulation is performed in two levels: (i) in the inner level we solve a series of *tightly coupled* three-dimensional problems in each subdomain using a semi-implicit parallel flow solver and (ii) in the outer level we solve explicitly *loosely coupled* problems by imposing proper velocity and pressure boundary conditions at the subdomain interfaces. We employ the continuous (*C*^{0}) Galerkin projection within patches and discontinuous-like Galerkin projection at the interfaces. The interface boundary conditions are based on the upwinding principles.

The potentially great advantage of the 2DD approach is not only in simulating very large-scale problems, such as the brain blood flow simulations that we consider in this paper, but also simulation of three-dimensional blood flow in the *entire* human arterial tree. In terms of computational efficiency, the scalability of the 2DD approach relies on the strong scaling and speed of a parallel solver employed at the inner level, and weak scaling at the outer level (Grinberg & Karniadakis 2008*b*). Depending on the problem size, our parallel solver, employed at the inner level, showed very good *strong* scaling on 1024–4096 processors of BlueGene and CRAY XT3, and hence we limit the number of processors assigned to each patch to 4096.

The 2DD method has been tested on the TeraGrid (Grinberg *et al*. 2007) and also on Ranger (http://www.tacc.utexas.edu/resources/hpcsystems). Our largest simulation to date (2.87 billions degrees of freedom) was performed on 18 576 processors and the solution was obtained at 1 s per time step. The computational domain, consisting of 147 arteries, was discretized into 1 244 295 tetrahedral elements and the solution was approximated with a sixth-order polynomial expansion. The simulation was performed to test the capabilities of the method and solver.

### (b) The RC method for pressure boundary conditions

The RC boundary condition for the pressure was developed as a *numerical method* to impose measured *in vivo* flow rates at multiple terminal branches. The method is valid for steady and unsteady flows, and it is based on the low-pass filter applied at each outlet to relate the flow rate and the pressure. The key idea of the method is in approximating the resistance (*R*) ratios at different branches by the inverse ratios of the corresponding flow rates (*Q*), i.e. *R*_{i}/*R*_{j}≈*Q*_{j}/*Q*_{i}, where *j* is the terminal branch ID. Then, the pressure at each time step is computed fromwhere *C*_{j} is a capacitance required to filter high-frequency oscillations, usually present in *Q*(*t*) due to numerical noise. In our simulations, the capacitance value was of order *C*_{j}≈0.2/*R*_{j} and the resistance *R*_{j} of order *O*(10)–*O*(1000). From the numerical and computational standpoints, the RC method is very efficient: it does not introduce additional constraints on the velocity system of equations and has very fast convergence, low computational cost and low volume of communications, which are very important in parallel computing.

### (c) PC-MRI measurements of flow rates

Cine phase-contrast imaging was performed with retrospective ECG cardiac gating on a 3T Trio system (Siemens Medical Systems, Erlanger, Germany). A slice perpendicular to the internal carotid and vertebral arteries was selected using a vessel scout as a reference, then a two-dimensional fast low-angle shot sequence (Frahm *et al*. 1986) was used in the slice of interest to form 25 images that represented sequential phases in a cardiac cycle. The following measurement parameters were used: flip angle 15°; TR=51.8 ms; TE=6.61 ms; slice thickness 5 mm; matrix size 512×512; 13.5 cm×13.5 cm field of view; and two signal averages. One of the acquired images is presented in figure 2. Velocity encoding of 60 cm s^{−1} allowed measurement without aliasing. The derived magnitude and weighted phase images were displayed in cine mode and qualitatively evaluated for flow. For each vessel, the region of interest (ROI) defined by cross-sectional area of its lumen was defined by thresholding. Volumetric flow rates (VFRs) in ml s^{−1} were then obtained from weighted phase images by integration of the velocities throughout the ROI. First, a zero-velocity baseline was estimated by calculating the average throughout of the cardiac cycle of the apparent velocity in a neighbouring static tissue. Second, VFR through the ROI was estimated by integrating the VFR through individual pixels (given by the product of the velocity in the pixel and the pixel area). Finally, VFR as a function of time waveform was calculated by repeating this procedure for each of the 25 frames. The derived waveforms for the left and right internal carotid and vertebral arteries, shown in figure 2, can be readily used in computational fluid dynamics (CFD) calculations as inflow boundary conditions.

## 3. Results

We performed two simulations of cranial blood flow using arterial geometry and inlet flow waves corresponding to two cases: (i) a healthy patient with complete CoW and (ii) a patient with hydrocephalus and incomplete CoW. Geometry of the two models was reconstructed from high-resolution computed tomography and MRI images using an in-house software package, a Matlab-based graphical user interface, that allows *accurate* detection of the arterial wall. Full automation of the geometry reconstruction process is not feasible, and in fact there is no commercial or free software that can guarantee automatic high-quality geometry reconstruction and meshing. It is crucial to control the segmentation process to avoid ‘digitally’ merged vessels, insertion of veins and bones as arterial walls. Our software was optimized for *interactive* processing of the images to manage the geometry reconstruction. The end product is a geometric model of the arterial wall saved in a PLOT3D or STL format and is compatible with different mesh generation softwares. We use Gridgen—a commercial mesh generator developed by Pointwise (www.pointwise.com), which also allows editing and construction of geometric databases. Typical time for geometry reconstruction is 2–12 hours, depending on the image quality and the user's experience. From the computational standpoint, the task can be performed on a standard desktop computer with 2–4 GB of RAM. The complete CoW arterial model, illustrated in figure 1, consists of 65 vessels with a diameter of approximately 4.5 mm at ICA, and approximately 1 mm at the terminal vessels. The computational domain, subdivided into four patches, was discretized with 459 250 tetrahedral spectral elements; in each element, the solution was approximated with polynomial order *P*=5. High-order spatial discretization required solution of four systems of equations, each with 180 026 000 unknowns (quadrature points), at every time step. The incomplete CoW, presented in figure 3, consists of 23 vessels, and was discretized with 130 962 tetrahedral spectral elements where the solution was approximated with polynomials of order *P*=4, which corresponds to 33 002 424 unknowns per each variable. In the first model, the right and left ICAs are interconnected with the basilar artery forming a complete circle, while in the second study we have a case where the left ICA is disconnected due to the lack of the communicating arteries. Hence, we concentrate on the flow interaction between the branches of the right ICA and the basilar artery. In both simulations, we have employed the rigid wall model. The inlet boundary conditions for the first simulation were obtained from velocity waveforms measured by Doppler ultrasound at the ICAs and the vertebral arteries, while the flow rates for the second simulation were obtained using PC-MRI measurements.

### (a) Three-dimensional simulations

We pursued two goals in the performed simulations: (i) to evaluate the feasibility of our new 2DD and RC boundary condition methods in real-life cases and (ii) to compare the flow patterns in the complete CoW model of a healthy subject and in the incomplete CoW model of a patient with hydrocephalus and to define goals for future investigations. From the numerical and computational standpoints, we have found that the 2DD method can be successfully applied for very large-scale flow simulations, and the RC boundary condition is an effective tool in imposing efficiently the prescribed flow rate ratios at multiple outlets. As was pointed out by Grinberg & Karniadakis (2008*a*), the accuracy provided by the RC model depends on the resistance values. These values are certainly patient specific, and additional care should be taken when we consider a pathological case. In our simulations, the resistance was of the order of 10^{10} Pa s m^{−3}, which is comparable with the data reported by Stergiopulos *et al*. (1992) and used by Alastruey *et al*. (2007) for one-dimensional modelling of cranial blood circulation. In the simulations of flow circulation in a patient with hydrocephalus, relatively high discrepancy between the prescribed flow rate ratios at posterior cerebral arteries (PCAs) and numerically computed ones (approx. 38%) was observed. To enhance the accuracy of the method, the resistance parameters must be increasing, which is consistent with the fact that patients with hydrocephalus have higher peripheral resistance than healthy subjects due to high intracranial pressure. Doubling the resistance values resulted in lower differences (17%) between that prescribed by the RC method flow rates and the computed ones.

Results of the flow simulation in the complete CoW are presented in figures 1 and 4. In figure 1, we plot flow rates and pressure drop at selected arteries *P*(*t*)−*P*_{ref}(*t*), where *P*_{ref}(*t*) is the pressure computed at inlet of the left ICA. In figure 4, the swirling flow is illustrated by instantaneous streamlines and vectors of the velocity field. The unsteadiness and phase shift in the flow waves imposed at inlets result in alternating flow direction in the communicating arteries. Such behaviour was also observed in one-dimensional simulations by Alastruey *et al*. (2007); however, the one-dimensional simulation cannot predict the secondary flows.

In simulation of the flow in the incomplete CoW, a relatively high pressure drop was observed through segment PCA-P1 (figure 3). This points to the need to perform a detailed study to understand its nature, and to answer the question whether the high Δ*P* is a consequence of the pathological processes and what is the role of the boundary condition modelling. High Δ*P* was also observed in other simulations performed using various settings for the outflow pressure boundary conditions. We concluded that the high Δ*P* is due to the phase shift in the flow rate waveforms imposed at the two inlets. This finding may affect the diagnosis of a vascular disorder in the cranial system. One way to validate our results is to perform a dye angiogram to compare the flow distribution in the cranial arteries with the results of the simulation. The dye angiogram requires X-ray scanning and is performed only when *specifically clinically* indicated due to two reasons: (i) a considerable amount of radiation and (ii) the invasiveness of putting catheters into the vessels of the leg and brain; this is a painful procedure and usually requires sedation and sometimes general anaesthesia, and carries a risk of stroke. A non-invasive and harmless technique to measure intracranial flow is Doppler ultrasound; however, in the case of our patient, this method is not applicable due to high skull-bone thickness. Currently, we are working on implementing high-resolution PC-MRI to obtain the three-dimensional velocity field in the major cranial arteries. This technique has the potential of replacing dye angiograms, but it requires variable velocity encoding to extract data from different vessels and different slices. Development of secondary flows in the incomplete CoW model was also observed. In figure 5, we plot the flow patterns close to the junction where the ICA bifurcates to the middle cerebral and anterior cerebral arteries, which is a typical region for developing aneurysms (Weir 2002).

### (b) One-dimensional simulations

The three-dimensional simulations provide a deep insight into flow patterns and wall shear stress and may be used for particle tracking instead of dye angiograms. However, sometimes the only information required is the direction of the flow, arterial pressure drops and flow distribution. Such requirements might be fulfilled by performing computationally inexpensive one-dimensional simulations. In order to compare the results of one- and three-dimensional simulations, we performed skeletonization of the model shown in figure 6; then, setting pressure at all outlets to constant and imposing the same flow rates at the two inlets, we performed one- and three-dimensional simulations. The one-dimensional simulation over eight cardiac cycles took only 20 min, while approximately 24 hours per cycle were required for the three-dimensional simulation on 256 cores of CRAY XT4. In figure 6, we compare the pressure drop between the terminal branches and the two inlets and the flow rates computed at terminal branches using the three- and one-dimensional models; the similarity in the waveforms is remarkable in most of the arterial domain. The discrepancies between the flow rates predicted by the one- and three-dimensional models usually appear at the vessels located farther from the inlets (outlets 1, 2, 11 and 12 in figure 6). The primary reasons for such deviation are: (i) the one-dimensional model does not take into account the angle at which a particular vessel bifurcates and (ii) the shape of the velocity profile at any point of the arterial tree is assumed to be the same, which introduces errors in modelling the wall friction in arteries with swirling or reversal flow.

## 4. Summary and outlook

In this paper, we discuss the modelling and computational complexity of intracranial arterial flow simulations and present results of two simulations in the MaN of the human brain. We overview two new numerical approaches developed specifically for this type of flow simulations, namely the 2DD and the RC outflow boundary condition. Simulations of flow in the complete CoW of a healthy subject and in the incomplete circle of a patient with hydrocephalus show development of strong secondary flows in the communicating arteries and ICA. The simulation in the incomplete circle predicts high pressure drop along the PCA-P1 artery, connecting the basilar with the internal carotid. Simulation of flow in the incomplete CoW requires setting the peripheral resistance to high values in order to ensure proper net blood flow to different arteries. The high pressure variation correlated with the phase shift at the inlet flow rates in the patient with hydrocephalus points to a new direction in diagnosing cranial vasculature disorders. The relatively accurate prediction of the flow distribution and pressure variation by one-dimensional modelling is encouraging and suggests that the three-dimensional modelling should be employed to calibrate the parameters assumed by the one-dimensional modelling. We envision that the one-dimensional models will be the default computational models to provide the first estimate of blood circulation and also in surgical planning. The role of the three-dimensional simulations is to complement the one-dimensional models and provide accurate information regarding three-dimensional flow patterns, wall shear stress, oscillatory shear index, low residence time zones, etc.

The ability to follow vascular lesions over time with less invasive means could be of clinical usage immediately. More sophisticated quantitative blood flow techniques may find near-term value in clinical medicine. We envision that CFD modelling of blood flow will replace many invasive and dangerous clinical tests.

While the present results are encouraging, more work is required to develop a truly multi-scale framework by coupling the three-level networks, namely MaN–MeN–MiN in a seamless integration. Ultimately, we aim at developing an efficient and scalable biomechanics gateway for analysing the biophysics of brain pathologies, such as hydrocephalus, in a systematic way on the emerging petaflops computer platforms of the next decade.

## Acknowledgments

This work was supported by the National Science Foundation (NSF) CI-TEAM grant OCI-0636336 and OCI-0845449; computations were performed at Texas Advanced Computing Center. E.C. acknowledges support under NSF award DUE-0734234. T.A. and J.R.M. thank the Webster family for supporting their work.

## Footnotes

One contribution of 15 to a Theme Issue ‘The virtual physiological human: tools and applications II’.

- © 2009 The Royal Society