Coupling one-dimensional arterial blood flow to three-dimensional tissue perfusion models for in silico trials of acute ischaemic stroke

An acute ischaemic stroke is due to the sudden blockage of an intracranial blood vessel by an embolized thrombus. In the context of setting up in silico trials for the treatment of acute ischaemic stroke, the effect of a stroke on perfusion and metabolism of brain tissue should be modelled to predict final infarcted brain tissue. This requires coupling of blood flow and tissue perfusion models. A one-dimensional intracranial blood flow model and a method to couple this to a brain tissue perfusion model for patient-specific simulations is presented. Image-based patient-specific data on the anatomy of the circle of Willis are combined with literature data and models for vessel anatomy not visible in the images, to create an extended model for each patient from the larger vessels down to the pial surface. The coupling between arterial blood flow and tissue perfusion occurs at the pial surface through the estimation of perfusion territories. The coupling method is able to accurately estimate perfusion territories. Finally, we argue that blood flow can be approximated as steady-state flow at the interface between arterial blood flow and tissue perfusion to reduce the cost of organ-scale simulations.


Introduction
An acute ischaemic stroke is the most common type of stroke. It is fatal to an estimated 3 million people globally every year and may lead to a dramatic loss of quality of life in survivors [1]. An acute ischaemic stroke occurs when blood supply to the brain is suddenly blocked by an embolized thrombus. This results in the acute loss of blood flow to downstream vessels. Consequently, brain tissue is starved of oxygen, resulting in necrosis, causing disability and ultimately death. Currently, only a few treatment options exist. Improving existing treatments and developing and approving new treatments is a difficult, lengthy and expensive process, prone to failures during clinical trials even after successful animal trials.
The INSIST project 1 [2,3] aims to develop computational models that will be used to simulate a clinical stroke trial on a virtual population. The (virtual) patient in this context is generally a set of parameters and other data that are used in a computational model. In silico trials (IST) can overcome some of the limitations of clinical trials, lower their cost, contribute to developing more efficient clinical trials or shorten the preclinical trajectories [4]. IST are closely related to computational biomedicine or personalized medicine. The goal is to predict the efficacy and efficiency of a treatment, drug or device. As such, they are similar to normal clinical trials where the testing happens on animals or humans. Computational modelling has become the standard in many industries from the design of air planes to electronics. Likewise, computational models in the medical field are slowly gaining momentum [5][6][7][8][9][10].
One of the goals of INSIST is to predict the location and volume of the infarcted brain tissue for stroke patients both for individual patients and at the population level. This requires modelling blood flow across length scales incorporating three orders of magnitude, from the large arteries via the arterioles and the pial surface vessels to the penetrating vessels and the microcirculation deep in the brain. Blood flow in large vessels is typically modelled using lumped parameter or one-dimensional (1D) blood flow models, whereas the microcirculation is typically modelled as a porous medium [11]. In addition, there is limited knowledge about the vessel anatomy between the large vasculature and the microcirculation. Data on the intermediate vessels are limited or only available as statistical scaling laws. Predicting infarcted volume requires thus the coupling of arterial blood flow models to tissue perfusion models. Inclusion of the collateral circulation is also important but is outside the scope of this paper. A tissue perfusion model is presented in another article in this same special issue of Interface Focus [12]. A 1D blood flow model has to be coupled to this tissue perfusion model, which presents several challenges.
Firstly, medical images often suffer from low resolution. Scans from stroke patients are often done quickly to get an assessment of the stroke severity. As a result, the scans typically only contain the circle of Willis (CoW) and a few of the major cerebral arteries. This presents a problem as this does not provide sufficient information for the coupling between arterial blood flow and tissue perfusion. However, there are high-resolution scans available from different individuals [13,14]. Unfortunately, this dataset, referred to as the BraVa dataset, does not contain other information such as, for example, pressure, heart rate and Young's modulus.
Secondly, information about the perfusion territories of the cerebral arteries can be determined but is limited [15]. The perfusion territory of a blood vessel is the region of tissue that receives blood from that vessel. Not much is known about the perfusion territories of the brain and data are again sparse. Coupling blood flow to tissue perfusion in a patientspecific manner requires the estimation of the perfusion territory for every perfusing vessel.
Here, a method is presented to couple blood flow in cerebral blood vessels to cerebral tissue perfusion. Patient-specific image-based vessel segmentations are combined with literature data and models for vessel anatomy not visible in the images to create an extended model for each patient. Using Murray's Law, we estimate perfusion territories on the pial surface using a coupling algorithm. Arterial blood flow is simulated with a 1D pulsatile blood flow model with three-element windkessel elements at the boundaries. We show that, at the pial surface, blood flow pulsatility is small and argue that blood flow can be approximated as steady-state flow.

Vasculature of virtual patients
Vessel centrelines are extracted from vessel segmentations of computed tomography (CT) angiographs from stroke patients entering the hospital [16]. 2 The centrelines contain the CoW and the side branches, figure 1a. The side branches are the posterior cerebral artery (PCA), anterior cerebral artery (ACA) and middle cerebral artery (MCA). This patient-specific vasculature is first extended with the larger vessels, starting from the heart, figure 1b. The smaller vessels to the cerebellum and brainstem are also added. The default parameters describing these arteries are taken from the literature [17,18] and are shown in table 1. A taper can be included by setting different proximal and distal radii but has not been considered in this paper.
Centrelines of high-resolution scans of cerebral arteries are available from 61 different individuals [13,14]. These centrelines are referred to as the BraVa dataset. The arteries in the datasets are labelled by six major cerebral arteries from which they emerged, i.e. left and right ACA, MCA and PCA, and the CoW, figure 1c. One segmentation of the BraVa set is randomly chosen and attached to the patient-specific vasculature. For each major cerebral artery in the patient vasculature, the equivalent vessel in the BraVa vasculature is used as an attachment point. The radii of the added vessels are scaled by the ratio of equivalent cerebral vessels, ACA, MCA, PCA, etc., between the BraVa dataset and the patient-specific vasculature. All radii are scaled to preserve any scaling law between the vessels.
These three sets of data are combined to create a patientspecific vasculature that starts at the heart and ends close to the pial surface. For each major cerebral artery in the patient dataset, the relevant vessels in the Brava dataset are identified and then attached to each other. On the other end of the patient anatomy, the large arteries to the heart are connected, as shown in figure 1d.

One-dimensional blood flow modelling
One-dimensional blood flow models assembled from multiple vessel segments [19][20][21][22][23][24][25] are able to model blood flow in sufficient detail without becoming too computationally expensive. These 1D models capture the axial flow component along the vessel, showing a good agreement with 3D blood flow simulations [26]. Blood flow is modelled as a viscous incompressible fluid. Assuming rotational symmetry around the centreline of the circular vessels and integrating the continuity equation along the azimuthal and the radial coordinates of a cylindrical coordinate system [19,22] leads to Following the same steps results in the averaged Navier-Stokes momentum equations The primary variables are the cross-sectional area of the vessel (A), the averaged axial velocity along the vessel ( v x ) and static pressure ( p) within each vessel as functions of the axial coordinate (x) and time (t). The corresponding model parameters are the momentum correction factor (α), the density (ρ) and the dynamic viscosity of blood (μ). A stroke is modelled as a complete blockage of flow and pressure by enforcing Neumann pressure boundary conditions at the end of the clotted vessel segment. In the future, the model will be extended to handle partial occlusions. In order to close the equation system formed by equations (2.1) and (2.2), a relationship between the pressure and the crosssectional area has to be defined. Wall deformations are described by a thin-walled axisymmetric elastic membrane approach so that : ð2:4Þ royalsocietypublishing.org/journal/rsfs Interface Focus 11: 20190125 Here, A 0 is the area of the lumen at p ref and p ref is a reference pressure set to the diastolic pressure. The β coefficient is determined by the vessel wall properties, which are constant in each segment, including the wall thickness (h), the Poisson's ratio (ν) and the Young's modulus (E) of the vessel wall (constants within vessel segments). The thickness of the vessel wall is calculated using an empirical function given by with r 0 the initial radius, a = 0.2802, b = −0.5053 mm −1 , c = 0.1324 and d = −0.01114 mm −1 [25]. The method used here to solve the resulting equations is the MacCormack method, a second-order finite difference method. The initial conditions are zero velocity and diastolic pressure everywhere. Vessels are linked together at the bifurcations by the method of characteristics for continuous propagation of characteristics, conservation of mass and total continuity of pressure (Bernoulli's equation) [19][20][21][22][23][24][25]. The resulting equations are solved with the Newton-Raphson method.
Every vessel has a minimum of three nodes and a maximum separation of 10 mm. The simulation is run for a small number of heart beats and iterated until the tolerance between iterations is below a threshold, defined as with p i being a vector of the pressure at all n nodes at all m timesteps during the ith cardiac cycle and p iÀ1 a vector of the pressure at the previous iteration i − 1. That is p i ¼ [ p 1 (t 1 ),p 1 (t 2 ) . . . p n (t m )]. The typical number of simulated cardiac cycles before convergence is below 10.
The heart provides the inlet boundary conditions to the model while a three-element windkessel model is used at each outlet. Note that the BraVa anatomies typically have of the order of 100 The large arteries, starting from the heart leading up to and including the patient-specific CoW. Smaller vessels to the brainstem and cerebellum are also included. Note that the vessels are not drawn to scale. (c) One of the centreline segmentations from the BraVa dataset (ID: BG18). Each vessel is smoothed by setting the radius at every point to the mean radius of that vessel. Vessels are colour coded according to which major cerebral vessel they belong. The CoW vessels from the BraVa segmentation are not used and are replaced by the patientspecific CoW from (a). (d ) The final patient-specific arterial network. Patient-specific data from (a) are extended with literature values from (b) and high-resolution data from (c). This is the network that is used in the 1D blood flow simulations. Note that the vessels are not drawn to scale.
royalsocietypublishing.org/journal/rsfs Interface Focus 11: 20190125 outlets. Volume flow rate at the inlet is given by an inlet function and taken from Boileau et al. [23]. The inlet function is scaled based on stroke volume and heart rate for each patient. The total resistance (R tot ) and compliance (C tot ) of the full vascular tree model are calculated by [27] with P sys the systolic pressure, P dia the diastolic pressure, SV the stroke volume, HR the heart rate and τ the aortic pressure decay time constant. The compliance of the patient network (C 1D ) is calculated from the sum of the vessel compliances given by where L is the vessel length. The network compliance is substantial for a large vessel and is, therefore, subtracted from the total compliance. The resistance of the patient network is assumed to be negligible for the large arteries compared to the total resistance. Known distributions of cardiac output into the main vessels are used to distribute the resistance and compliance to each outlet. That is, 65% goes to the lower body through the thoracic aorta, 5% goes to the arms through the left and right brachial arteries and the rest to the other outlets [28]. The resistance and compliance are distributed to the outlets, R T and C, based on the fraction of the radius cubed, e.g. Murray's Law [29]. The cardiac output distribution is assumed to be constant. The total resistance and compliance for an outlet are given by where i sums over all outlets connected to a particular body part and CO f is the cardiac output fraction of that body part. The total resistance of the attached cerebral arteries is subtracted to account for the resistance added by the attachment of the BraVa vessels to the patient network. The resistance of the added vasculature is calculated by iteratively adding resistances depending on whether vessels run in parallel or serial. The two resistance windkessel parameters, R 1 and R 2 , of each outlet are calculated by and The choice for R 1 leads to the minimization of unnatural wave reflections at the outlets [30]. The outlets of the 1D blood flow model provide blood flow to the perfusion territories of the brain. The default parameters of the model are listed in table 2.
The extension process proceeds along the following steps: (1) identify identical vessels in both datasets by their label; (2) extract the smaller upstream vessels starting from the major cerebral artery; (3) calculate the mean radius for identical vessels in both datasets; (4) scale the radii of the upstream vessels identified previously accordingly; (5) attach the upstream vessels to the relevant end in the patient vasculature; (6) calculate the resistance of the appended vessels; (7) subtract the resistance and distribute the remainder to the appended ends.
To characterize the strength of the pulsatility of blood flow, pulsatility indices are used. The pressure pulsatility index is given by with p max the maximum pressure, p min the minimum pressure and p mean the mean pressure over a single heartbeat. The velocity pulsatility index is given by with v max the maximum velocity, v min the minimum velocity and v mean the mean velocity.

Estimating perfusion territories
Coupling blood flow to tissue perfusion in a patient-specific manner requires the estimation of the perfusion territory for every perfusing vessel, so for every outlet of the BraVa sets (i.e. left and right ACA, MCA and PCA), the cerebellum and the brainstem. Using Murray's Law, perfusion territories on the pial surface are estimated using a coupling algorithm. The radius of the outlet is assumed to be related to the perfusion of tissue by Murray's Law given by with r 3 i the cubic radius at outlet i and Q i the volume flow rate at outlet i.
We estimate the fraction of the pial surface that belongs to each outlet i of the patient vasculature at the pial surface with the use of Murray's Law. This fraction is calculated as with Q i the volume flow rate at the outlet and r 3 i the cubic radius at the outlet. The proportionality constants from equation (2.16) drop out in equation (2.17) because they are set to be the same at all outlets. The volume flow rate per area, or flow velocity, is uniform over the pial surface under these assumptions. To what extent this is also true for a real brain is not clear.
The pial surface is represented as a uniform triangulated mesh, figure 2a, where each triangle represents an equal sized area on the mesh. The mesh was presented in the study of Garcia-Gonzalez et al. [31] where it was used for finite-element computations. The triangles are grouped based on distance, calculated with Dijkstra's algorithm, using a coupling algorithm, figure 2b and c. The number of surface elements (n i ) that belong to each outlet is calculated by with N Tot the number of triangles of the uniform triangulated mesh.
The coupling algorithm ensures that the size of the perfusion territories is scaled based on the cubed radius of vessel that perfuses it, i.e. Murray's Law. After the coupling, any unassigned surface elements are assigned to the nearest cluster. The surface elements are assigned to the nearest outlet based on their separation on the surface. This is an iterative algorithm to account for the mismatch between the pial surface and the outlet. The outlet position is updated by taking the centre point of the assigned surface elements as a new guess. The new centre point is determined by minimizing the maximum separation between points.
Data from vessel-encoded arterial spin labelling (VEASL) is used to confine the sampling of the mesh to regions that belong to the same major cerebral vessel [32][33][34]. This is done to achieve a better match between the estimated area and area of the surface mesh. The pial surface mesh does not contain the complex folding of the brain nor does it contain well separated frontal and temporal lobes. By using VEASL data, the mismatch between the expected royalsocietypublishing.org/journal/rsfs Interface Focus 11: 20190125 area and the area of the surface mesh can be limited to the major territories of the brain. The major territories of the brain are the left and right ACA, MCA and PCA, the cerebellum and the brainstem. Vessels perfusing the pial surface are confined to the equivalent regions in the VEASL data. The coupling algorithm tries to find the sets for which the sum of distances within the set of surface elements is minimal. Note that the cluster sizes are calculated beforehand using equation (2.18). The root of the perfusion territories can be thought of as the mean position of the surface elements in that set. The surface is transformed into a graph and Dijkstra's distance algorithm is used, figure 2b. Euclidean distance can also be used if the surface is relatively flat. The pseudo code is given as follows: (1) Determine the number of surface elements per outlet with equation (2.18). (2) Project the outlet to the pial surface as an initial guess for the root. The mean flow velocity at outlets of the blood flow model is distributed evenly across the relevant perfusion territory. Each triangle of the mesh receives a volume flow rate based on its fractional area such that the flow velocity, i.e. volume flow rate per area, is uniform within each perfusion territory and within the major region set by the VEASL data. The result of the coupling is shown in figure 4, error estimation in figure 5 and blood flow at the pial surface in figure 6.

Results and discussion
Extending patient-specific data on the anatomy of the CoW with high-resolution data from the BraVa sets and literature can solve some of the problems of missing vessels. We include the arteries starting at the heart such that we can capture the redirection of flow in the event of a stroke. This is often excluded from many models that look at stroke. For large vessel occlusions, it is likely that the redistribution of blood can have significant effects on perfusion. A more detailed study on this is in progress and will be reported elsewhere.
The extended patient vasculature, figure 1d, should contain sufficient detail to capture all relevant effects such as the redistribution of flow after a stroke and make it possible to model volume flow rates, pressure and other variables. Once the patient anatomy is extended to the smaller cerebral arteries, we can begin to map the flow in the arteries to the brain. Modelling the entire cerebral vasculature down to the microcirculation is computationally unfeasible and probably also not required for our purpose. Being able to switch to a more efficient model is, therefore, crucial. The 1D blood flow model used in this paper is commonly used, well studied and validated [16,23,25]. The outlets of the 1D blood flow model provide blood flow to the perfusion territories of the brain.
Statistically accurate vessel networks can be generated between the outlets of the BraVa sets and the coupling points at the pial surface (figure 1c) using constrained constructive optimization [35][36][37]. However, one of the underlying assumptions of such algorithms is that the radii at the terminal points are equal and their positions are randomly chosen within a target volume. Since the assumption here is that there is no seepage through the vessel wall, the volume flow rates at the outlets would all be equal for a symmetric bifurcating tree. Simple scaling laws, such as Murray's Law and a length-toradius ratio, can be used to generate bifurcating trees, as the 1D blood flow model only depends on the length and the radius and not the position of the vessels.
Pulsatile blood flow was simulated in the anatomy of figure 1d, using parameters from table 2. The simulations were performed for seven heartbeats, after which the solution converged. Time-dependent volume flow rate and pressure along the anatomy were obtained, and from that both velocity and pressure pulsatility indices (PI) were computed.
The pressure and velocity pulsatility decrease as a function of distance from the heart for most vessels. Pulse pressure amplification is observed at the brachial arteries. Figure 3a shows the volume flow rate, pressure and radius of three vessels at different sizes, one large, one intermediate and one small artery. The chosen vessels were the right common carotid artery, the right MCA and one outlet vessel belonging to the right middle cerebral vasculature. The figure shows that the volume flow rate profile only changes in magnitude as the flow is split at every bifurcation. The pressure profile changes as the higher frequencies are damped. The vessel radius changes depending on the pressure and the properties of the vessel. The radii of the small cerebral arteries do not change much. The cerebral vessels are much stiffer than the large arteries and the pressure is lower further away from the heart, figure 3b [38]. The cerebral vessels beyond the CoW are similar in size; their Young's modulus is assumed to be similar and is set to the same value as their main branch of 1.6 MPa. Figure 3c shows the pressure PI in the cerebral vasculature. Figure 3d shows the pressure and velocity PIs as a function of distance from the heart. Measurements of the velocity pulsatility index of the right MCA are possible through transcranial Doppler ultrasound. The PI measured is often the velocity PI also known as Gosling's pulsatility index [39]. PI is thought to reflect the downstream vascular resistance and changes in PI are associated with various diseases [40,41]. Table 3 shows the velocity pulsatility indices for the MCA. Comparing the values found to those in the literature shows that the model overestimates the velocity PI [40,[42][43][44]. Stiffness of the vessels plays an important role in the transmission of pressure and velocity pulsatility through the vasculature [44,45]. Obtaining patient-specific Young's moduli for every vessel is difficult and can explain some of the overestimation of the velocity PI in the model. The outlet resistance and choice of BraVa segmentation affect the pulsatility indices more than the Young's moduli. Pressure pulsatility index decreases more than the velocity pulsatility index as a function of distance from the heart. This is likely due to the dampening seen in the pressure profiles, figure 3a. This suggests that the pressure PI is a better measure of the downstream resistance than the velocity PI.
The pressure and velocity pulsatility decrease as the radius of the vessel decreases, figure 3c and d. High frequencies are effectively decreased in magnitude. The pressure pulsatility index shows that there are two types of behaviour. The first type is the pulse pressure amplification in the brachial arteries. The second type is the slow decrease in pulsatility index as distance from the heart increases. Velocity pulsatility seems to reach a minimum at around one. Tables 3 and 4  royalsocietypublishing.org/journal/rsfs Interface Focus 11: 20190125 used to personalize the model. The configuration of the CoW and the stroke itself will change the outcome such that we cannot compare the model output to values from the literature for healthy patients. Note that the model is not calibrated nor fitted to match the reference values. The pulsatility indices when compared with the reference value do seem to be too high but, without either calibration or sufficient patient data, this is not surprising. The configuration of the CoW plays a major role in cerebral blood flow, especially during a stroke [16,46]. Including more patientspecific parameters, such as length, radii, elasticity, etc., will affect the model results. However, these are often difficult or infeasible to obtain for stroke patients. Getting correlated high-quality data for stroke patients is difficult. This is something that we hope to address in future work. The coupling algorithm is able to consistently map the outlets of the blood flow model to the pial surface. The coupling method works with any pial surface mesh and a higher-resolution mesh will also increase the resolution of the perfusion territories. The choice of BraVa subject does not affect the blood flow variables before the CoW with the coefficient of variation not exceeding 0.005 for any vessel in table 1. However, there is variation in the location and size of the perfusion territories. This likely reflects the individuality of the cerebral vasculature. Computation times for the 1D blood flow simulation are of the order of 10 min while the coupling algorithm takes significantly more time and is of the order of 60 min. Computations are done on a standard desktop computer with an Intel core i7-7700 k running at 4.2 GHz with 16 GB RAM. Figure 4a shows the mapping of the pial surface based on the major artery from which the vessels originate. To date, only the perfusion territories for the major cerebral arteries are well established [15,47]. The major cerebral perfusion territories are enforced with the use of VEASL data. This is done to try to limit the error between the expected area and the area on the surface mesh. The surface mesh does not capture all the complex folding of the brain and a mismatch is expected. This lack of resolution is a serious limitation of the coupling algorithm. A downside of this is that the fractional areas, equation (2.17), are only strictly valid within each region, figure 5.
The coupling algorithm essentially maps the pial surface to the nearest vessel. The main benefit from the coupling algorithm is that it ensures that the perfusion territories are sized correctly, figure 5. Shown in figure 5 are ratios between various variables within regions of the brain. Every outlet is compared to every other outlet in a specific region of the brain. The ratios provide a way to compare the results of the model to the model assumptions as a form of verification. Figure 5a shows the relationship between the cubic radius and the area on the pial surface of the vessel outlets. Figure 5b shows the relationship between cubic radius and mean volume flow rates, see equation (2.16). Figure 5c shows the relation between the number of triangles and area on the surface mesh. The ratios display a linear relationship; any deviation from the linear relationship represents a deviation from the model assumptions. The deviations come from the fact that the resistance of the vasculature affects the flow and pressure at the outlet, and the surface segmentation is divided into a discrete number of triangles (numerical error).
The coupling method is graph-based and works with any geometry and distance metric. The coupling method uses equation (2.17) to estimate the fractional area of each outlet without the need to calculate flow rates. Once the outlets are mapped to the surface, the effect of different diseases, such as a stroke or stenosis, on blood flow can be modelled.
The assumption is that the underlying vasculature does not change due to the disease. Different scaling laws can be applied by changing the calculation of fractional area, equation (2.17). Using a different scaling law affects the size and location of the perfusion territories. It would be interesting to see the effect of different scaling laws. Murray's Law is royalsocietypublishing.org/journal/rsfs Interface Focus 11: 20190125 used in this paper as it seems to be a good fit for the cerebral vasculature [48]. A simple alternative to the coupling algorithm would be to map the surface to the nearest outlet. However, this can cause problems when some outlets are mapped to large areas of the surface and others are not mapped at all. The effect of differently sized radii at the ends of the vessels cannot be ignored. Another way to estimate perfusion territories is based on calculating streamlines [49]. However, this method is computationally expensive and can run into issues with large complicated geometries.
Changes in volume flow rate become negligible on the time scale of tissue oxygenation and cell death. We, therefore, argue that we can approximate the volume flow rates at the pial surface as steady-state flow. The mean volume flow rate over the region does not depend on details of the vasculature due to conservation of mass. The steady state only changes due to a major event in the vasculature such as a stroke. We, therefore, map the mean volume flow rate directly from the vessel outlet to its perfusion territory. One could go further and replace the entire blood flow model with a steady-state resistance model.
The mapping of the pial surface to the outlets of the 1D blood flow model is shown in figure 4b. The volume flow rate at the outlets of the 1D blood flow model is distributed evenly across its perfusion territory. Figure 6a shows the mean blood flow velocity at the pial surface and figure 6b shows the mean volume blood flow rate at the pial surface. There are some differences in flow rate at the pial surface but all flow rates are within one order of magnitude. These differences are largely due to deviations from the assumed proportionality between radius and volume flow rate, equation (2.16) and figure 5b. The flow rates at the surface depend not just on the radius at the outlet but also on the resistances up to the outlets. The surface mesh also does not capture all the complexity of the human brain and affects royalsocietypublishing.org/journal/rsfs Interface Focus 11: 20190125 the distribution of the flow across the pial surface, figure 5a and c. There is a trade-off between increasing the number of surface triangles to improve accuracy and computation time. Improving the quality of the surface mesh and improving the surface mapping as well as validation of the mapping will be addressed in future work. An extensive validation based on the MR CLEAN database is in progress. Infarcted volume prediction based on the model presented here coupled to a porous model shows promising agreement with CT scans [12].
A clot can be modelled by introducing a reflective boundary condition at the boundary of the flow. This is achieved by enforcing a zero pressure gradient, resulting in zero flow through the clot. The redistribution of flow after a stroke is included in the model. The flow to the affected region is redistributed to the other outlets. The brain receives a small portion of this flow; most goes to the other parts of the body. The brain receives only about 15% of cardiac output. A first approximation would, therefore, be that the brain receives 15% of the flow to the occluded region. Blood flow to the rest of the brain is not expected to be affected significantly. Figure 6c shows the effect of a clot in the right MCA on the velocity of flow at the pial surface. Figure 6d shows the effect of the same clot on the volume flow rate to the pial surface. The figures show that downstream of the clot there is no flow, leading to regions on the pial surface that lack flow. These are the areas where an infarct will start to form.
The framework presented here can be used to model blood flow to organs such as the brain. In addition, the  Figure 6. (a) Mean blood flow velocity at the pial surface for a healthy patient. The mean volume flow rate divided by the cluster area for each vessel outlet to the pial surface is shown. All values are within one order of magnitude of each other. (b) Volume flow rate on the mesh for a healthy patient. The volume flow rate at the outlet of the vessel is mapped to its cluster on the mesh. There are large differences between clusters, as the outlet radii can differ significantly. (c) Flow velocity at the pial surface for a patient with a stroke. There is a clot in the right MCA. Pressure after the clot is equal to that of the venous system and flow is zero. (d ) Volume flow rate at the pial surface for a patient with a stroke. Any vessel originating from the right MCA has zero flow and their perfusion territories receive no blood. coupling algorithm can be used to estimate perfusion territories for any organ as long as the appropriate distance metric is used. It is important to capture enough of the vasculature to be able to accurately predict pressure and flow rates. Obtaining correlated parameters as well as patientspecific segmentations is a major challenge for creating patient-specific simulations.

Conclusion
In silico trials can dramatically improve the process of development of medical devices, drugs and treatments. IST can potentially reduce the cost and duration of running a trial by reducing the need for patients and speeding up testing. The INSIST consortium [50] aims to create an IST for treatment of acute ischaemic stroke. The trial will include models of the patient population, arterial blood flow, tissue perfusion, metabolism, thrombosis, thrombolysis and thrombectomy. Here, the arterial blood flow model and a method for estimating perfusion territories are presented. Imagebased patient-specific data on the anatomy of the CoW are combined with literature data and models for vessel anatomy not visible in the images to create an extended model for each patient. An acute ischaemic stroke can be modelled as the blockage of flow at any point in the vasculature. Modelling stroke and predicting infarct volume requires the coupling of multiple models on different scales. A method to couple patient-specific blood flow models to a tissue perfusion model is presented. Blood flow to the pial surface can be approximated as steady if we accept a small error due to the loss of the time-dependent behaviour. However, additional work remains to be done on the incorporation of feedback between the blood flow model and the tissue perfusion model.