Abstract
There is a growing need for patientspecific and holistic modelling of the heart to support comprehensive disease assessment and intervention planning as well as prediction of therapeutic outcomes. We propose a patientspecific model of the whole human heart, which integrates morphology, dynamics and haemodynamic parameters at the organ level. The modelled cardiac structures are robustly estimated from fourdimensional cardiac computed tomography (CT), including all four chambers and valves as well as the ascending aorta and pulmonary artery. The patientspecific geometry serves as an input to a threedimensional Navier–Stokes solver that derives realistic haemodynamics, constrained by the local anatomy, along the entire heart cycle. We evaluated our framework with various heart pathologies and the results correlate with relevant literature reports.
1. Introduction
The strong anatomical, functional and haemodynamic interdependency of various cardiovascular structures are nowadays greatly appreciated. An increased holistic view of the cardiac system demanded by clinicians is especially valuable in the context of heart disease and associated complex dysfunctions. This is in perfect accordance with the tremendous scientific effort worldwide, such as the Virtual Physiological Human [1] project geared towards multiscale physiological modelling and simulation, which will promote personalized, preventive and predictive healthcare.
In recent years, computational fluid dynamics (CFD)based techniques have been used with varying degrees of success for simulating blood flow in the heart [2–7]. Most of these techniques employ a moving boundary setup, where the motion of the walls is prescribed by extracting the dynamics from the fourdimensional imaging data. Recently, there has been work reported on fully coupled fluid–structure interaction (FSI) techniques [8], which, in addition to CFD simulations, also solve the coupled solid mechanics problem for the wall deformation. One of the main challenges of FSI techniques is the characterization of the advanced material property model to accurately reflect the properties of the cardiac tissues.
For a realistic patientspecific haemodynamic simulation, one not only needs an accurate modelling of the computational mechanics (solid and/or fluid), but also the comprehensive anatomical model of a patient's heart, containing the four chambers, valves and main vasculature. In order to address these issues, we propose an integrated model of the entire heart, including anatomical, dynamical and haemodynamic representations (figure 1). Part of this work has been reported in several of our conference publications (e.g. [9]). In this paper, the novel contributions include: complete modelling of the right heart chambers, valves and vessels; an extended account of our numerical haemodynamics computation setup; extended description and analysis of our new haemodynamics simulations, including vortical dynamics and the influence of the atria.
This paper is organized as follows. In §§2 and 3, we describe the machinelearning algorithms used for robust estimation of the patientspecific parameters of morphology and function from fourdimensional cardiac computed tomography (CT) data. Based on the dynamic model of the complete heart, CFD simulations (in a moving boundary framework) are employed to simulate the cardiac blood flow. The details of the CFD modelling are discussed in §4. In §5, we present results that we obtained by applying our comprehensive approach on multiple fourdimensional CT datasets. We conclude the paper in §6, by revisiting the main ideas and identifying directions for future work.
2. Comprehensive heart model
The proposed comprehensive heart model includes the four cardiac chambers, valves, great vessels and captures anatomical, dynamical and pathological variations. To facilitate effective estimation of patientspecific parameters, the inherent complexity of the complete heart model is represented on two abstraction levels: landmark model and surface model. The landmark model represents the key location of critical anatomical locations as well as their synchronous nonlinear motion during the cardiac cycle. Each landmark L_{n} is parameterized as a T timestep trajectory in the threedimensional Euclidean space: L_{n} = {l_{1}, l_{2}, … , l_{T}}, n ∈ 1 … 43, l_{i} ∈ .
In order to capture the key anatomical and dynamical variations, we include the 43 landmarks, defined as follows: left ventricle apex, right ventricle apex, four landmarks for approximate location of the pulmonary veins, two landmarks for the venae cavae, two for the coronary ostia, seven landmarks for the mitral valve annulus and leaflets, 11 landmarks for the aortic root and leaflets, nine landmarks for the pulmonary valve and six for the tricuspid valve.
The surface model is constrained by the previously defined landmarks, and those together build the proposed comprehensive model. Each anatomical structure is represented separately by a surface mesh spanned along two physiological directions: M_{q}(L) = {ν_{1}^{q}, ν_{2}^{q}, … , ν_{K}^{q}}, q ∈ 1, … ,18, ν_{i} ∈ , where ν_{i} are the vertices, and K^{q} is the total number of vertices of a certain surface mesh q. A rectangular grid with n × m vertices is represented by (n − 1) × (m − 1) × 2 triangular faces. The complete model of the heart includes the following surfaces: left ventricle endocardium and epicardium, right ventricle endocardium, left atrium with corresponding pulmonary veins, right atrium with attached vena cava and the aortic, mitral, tricuspid and pulmonary valves.
To enable the construction of statistical shape models, point correspondence is maintained between models from different cardiac phases and across patients. The consistent representation is achieved through a set of physiologically defined sampling schemes, parametrized by the anatomical landmarks. For instance, for the aortic root, a pseudoparallel slice method is used, which equally distributes a sequence of cutting planes along the corresponding centreline (figure 2). To account for the bending of the aortic root between the commissures and hinge levels, specific to root dilation pathologies, at each location, the cuttingplane normal is aligned with the centreline tangent. The resulting twodimensional contours can be uniformly resampled to establish the desired point correspondence.
3. Cardiac anatomical and dynamical computation
The patientspecific parameters of the model are hierarchically estimated from fourdimensional cardiac CT images by employing a robust learningbased framework [10,11]. The a posteriori probability p(L, MI) of the models L, M given the image data I is incrementally detected within the marginal space learning (MSL) [11] framework. Detectors are successively trained using the probabilistic boosting tree (PBT) [12], with Haar and steerable features, and subsequently applied to estimate the anatomical landmarks and structures from cardiac CT volumes.
Given a sequence I of T volumes, the task is to find the model parameters L, M with maximum posterior probability: argmax_{L,M}p(L, MI) = argmax_{L,M}p(L_{1}, … , L_{41}, M_{1}, … , M_{18}I(0), … , I(t)). To solve the equation, we formulate the model estimation problem as a classification problem, for which a learned detector D models the probability p(L, MI), and scans exhaustively all parameter hypotheses to find the most plausible values for L and M. However, a direct estimation is computationally unfeasible, as the number of hypotheses H to be tested increases exponentially with the dimensionality of the model parameters. Assuming a simple model with only 10 parameters, each discretized to 10 values, the total numbers of hypotheses H is 10^{10}. To overcome this limitation, the MSL framework breaks the original parameter space Σ into subsets of marginal space with increased dimensionality Σ_{1} ⊂ Σ_{2} ⊂ … ⊂ Σ_{n} = Σ, where dim (Σ_{1}) < dim (Σ) and dim (Σ_{k}) − dim (Σ_{k}_{−1}) is small. A search in Σ_{1} with a detector D_{1} learned in this marginal space, finds a subspace C_{1} ⊂ Σ_{1}, which contains only the most probable parameter values, and discards the rest of the space, such that C_{1} < H_{1}. Another stage of training and testing is performed in the extended space C_{e}_{1} = C_{1} × H_{2} ⊂ Σ_{2}, to obtain a restricted marginal space C_{2} ⊂ Σ_{2}. The procedure ends when the final dimensionality of Σ is reached. Instead of using a single detector D, we train detectors for each marginal space and estimate the model parameters by gradually increasing dimensionality. After each detection stage, only a limited number of highlyprobable candidate hypotheses is maintained and considered as input for the following stage. The comprehensive model of the heart (containing the four chambers and the valves) thus obtained is illustrated in figure 3. For more details please refer to Ionasec et al. [13].
The performance of the model estimation approach was evaluated on a total of 1013 cardiac CT volumes from 206 patients affected by various cardiovascular diseases. The ground truth for training and testing was obtained through a manual annotation process, which was guided by clinical experts. The model estimation accuracy of the complete heart was measured using the pointtomesh distance to compare the ground truth and automatically detected models. An average accuracy of 1.36 ± 0.93 mm was obtained for the valves [13] and 1.13 ± 1.57 mm for the chambers [11], at a computation time of 8.8 s per each volume. It is important to note the robustness of the valve estimation method, with a median of 1.30 mm and 80percentile of 1.53 mm. Moreover, the precision of the patientspecific valve model estimation lies within 90 per cent of the interuser confidence interval.
4. Cardiac haemodynamics computations
4.1. Blood as a fluid continuum
The haemodynamics computations presented here use a classical continuum model for the blood, and are performed using the framework presented in Mihalef et al. [9,14]. The incompressible Navier–Stokes equations with viscous terms (4.1), which are the standard continuum mechanics model for fluid flow, are solved using direct numerical simulation in a level set formulation [15]: 4.1
Navier–Stokes are partial differential equations describing the momentum and mass conservation for fluid flow, depending on the velocity u and pressure p of the fluid, as well as on the fluid density ρ and dynamic viscosity μ. The blood density and dynamic viscosity are set to generic mean values across healthy individuals, namely ρ = 1.05 g cm^{−3} and μ = 4 m Pa · s. The equations are solved here using numerical discretization on a uniform grid, employing both finite difference and finitevolume techniques. In particular, we use the fractional step combined with an approximate projection method for the pressure proposed by Bell et al. [16].
We model the blood as a Newtonian liquid. Previous numerical studies (e.g. [17]) have found that, in larger arteries, the nonNewtonian behaviour is not important during most of the cardiac cycle. The nonNewtonian importance factor (defined as the ratio of the nonNewtonian and the Newtonian viscosity of blood) becomes significant during a 15–20% subperiod of the cardiac cycle, achieving a peak during a speed deceleration period when velocity is close to zero. Blood dynamics in the heart is mostly governed by high velocities—or, rather, high shear rates. This, combined with the fact that blood's rheology is close to that of a Bingham plastic liquid (a material that behaves as a rigid body at low stresses but flows as a viscous fluid at high stress), supports the qualitative conclusion that during the heart cycle, blood's dynamics is predominantly Newtonian. The probable exceptions are: timewise the diastasis, and spacewise the apical and the jet stagnation regions. A quantitative analysis of such phenomena would be best performed using enhanced models of the ventricle, including at the least the papillary muscles, and is beyond the scope of this paper.
4.2. Computational fluid dynamics in the level set framework
Computing fluid dynamics in moving generic geometries and/or involving multiphase flow constitutes a challenge, especially when it comes to devising computational methods that are simple and robust, as well as accurate. The level set methods [18] have achieved success in the recent years in dealing with such complex computations, being able to capture the intricate dynamics of interfaces between different materials (for example, water and air, or a complex deforming object and a liquid). The Navier–Stokes equations are cast into a form that takes into account such a description, for example, for two fluids with different density and viscosity, we get: 4.2In equation (4.2), H is the Heaviside function, used to create a numerically sharp distinction between the first fluid, characterized by positive values of the level set φ, and the second fluid characterized by negative ones. The level set formulation, especially when implemented using advanced methods like the ghost fluid method [19], has several advantages over classical formulation, for example, ease of computations of free surface flow and deforming materials interacting with fluids, simple implementations for tensor extrapolation technologies, and also simple formulae for various geometric parameters, like normal fields () or mean curvature fields (for example, κ = ∇ φ for level set functions of unit gradient).
4.3. Levelset modelling of the heart
Another important advantage of describing the heart dynamics in the levelset framework is the possibility of automatically dealing with moving/ heavily deforming walls without the need to generate a new computational grid at every several time steps, when the previous grid becomes too skewed to handle robustly the numerical computations, as is the case in finiteelement methods (FEM). This also allows for automatic integration of patientspecific cardiac meshes into the blood flow simulation engine, and offers a framework for clinical usage of cardiac models for obtaining patientspecific haemodynamics.
In this work, we model the heart walls/blood interface using a level set. More precisely, we start from the heart mesh sequence (comprising 10 frames), obtained as described in §§2 and 3, and interpolated using splines to derive the mesh at the given simulation time. This mesh is embedded in the computational domain with the help of the levelset function φ, defined as φ(x) = dist(x, mesh) − dx, where dx is the grid spacing. The heart/blood interface as used in the code is defined as the zero level of the levelset function, effectively ‘thickening’ the original triangle mesh by dx on each side (figure 4).
The interface location is used to impose noslip boundary conditions to the fluid region (see next section). The mesh velocity at each time step is easily computed by temporal interpolation from the mesh positions at adjacent time steps, then extrapolated to the interfacial nodes using extrapolation kernels. Note that imposing the boundary conditions in this manner effectively enforces oneway transfer of momentum from the solid heart mesh to the fluid. This approach essentially models the heart as a pump pushing against the domain walls, which can approximate the resistance of the circulatory system.
4.4. Algorithmic and numerical details
This section describes in more detail the numerical solution of the Navier–Stokes equations (4.2). We solve all stages of the cardiac cycle, including the isovolumic phases (IP). Note that recently Sengupta et al. [20] used highresolution Doppler to reveal that the IP are not periods of stasis but rather phases with dynamic changes in the intracavitary flow.
Our algorithm starts at a given time step n from the level set, velocity and pressure information at the previous time step φ^{n}^{−1}, u^{n}^{−1}, p^{n}^{−1}, and computes φ^{n}, u^{n}, p^{n} following a fractional step projection method:

— convective updates for φ and u;

— semiimplicit update for u to take into account the viscous force contribution;

— pressure update by solving the Poisson equation with Neumann boundary conditions; and

— new velocity update.
In step 1, we update first the levelset values, as described in §4.3, by using the mesh location at the given time step. In an intermediate geometric step, we compute the connected components defined by the new level set. This is used subsequently for robust connected componentbycomponent inversion of the implicit viscous and pressure Poisson linear systems. Then the convective force terms are computed using thirdorder accurate essentially nonoscillatory (ENO) techniques [21]. In the second step, the velocity is updated to u* by using the secondorder semiimplicit splitting proposed by Li et al. [22], and by inverting the system using an efficient multigrid preconditioned conjugate gradient solver. The same routine is used in step 3 to solve the pressure Poisson equation in each of the connected components of the discretization domain. Before the system inversion, we overwrite the velocity in the solid regions using the solid velocity. In step 4, we update the density using the new level set and compute the velocity update u^{n} = u* − ∇p^{n}/ρ^{n} in the liquid, or u^{n} = u^{solid} in the solid.
The interface location is used to impose noslip boundary conditions to the fluid region in, namely u_{fluid} = u_{solid} for the convective and the viscous component of the Navier–Stokes momentum equation, and ∂p/∂n = (u* − u_{solid}) for the pressure Poisson equation (here is the normal vector field, computed from the level set as described in §4.2). The global accuracy of the method is secondorder in the interior of the domain, and this drops to first order at the boundaries. For further details concerning the formulation and the solution of this approach, please refer to Mihalef et al. [14,23], and Sussman [24].
4.4.1. More about boundary conditions
The human heart functions as a part of the cardiac system, thus being influenced by the loading pressures from the venous and arterial system. Our simulation framework solves the Navier–Stokes equations with prescribed (wall) boundary motion, thus ensuring that we do not need to compute the pressure field but rather its gradient, as the pressure is a relative rather than absolute variable. In other words, the first Navier–Stokes equation in equation (4.2) (momentum conservation) is invariant to gauge shifts of the pressure, as long as its gradient does not change. Mathematically, the pressure Poisson equation with Neumann boundary condition has a oneparameter family of solutions, which can be fixed once a Dirichlet boundary condition is introduced. Therefore, when solving the pressure Poisson equation, we fix the solution by choosing a base pressure (equal to zero) at a point outside the aorta, and this acts in a similar way with imposing p = 0 at the aortic outflow surface.
As a note, we can choose any value for the base pressure, for example, using the values from a generic physiological curve, and this does not influence the haemodynamics in the prescribed motion framework presented here. This would however influence the dynamics if one would solve FSI problems, where deformable tissues would respond to pressure boundary value changes, or in the case of regurgitations. For such cases, the coupling of our model with a onedimensional systemic vascular model would be necessary, and we are working towards such integration.
5. Results
Using the framework described in §4, we performed a series of blood flow simulations using the patientspecific heart model. By measuring geometrically the stroke volume in the left and right ventricle we obtained values that agreed within a maximum difference of 5 per cent. The left and right sides were not connected through the systemic circulation, owing to the lack of geometric and material data for a full vessel tree. As a consequence, we performed the simulations separately for each of the sides of the heart, which is computationally more efficient, and also is not influenced by the slightly different stroke volumes in the left and right sides.
Snapshots of our simulation results are visualized in figure 1. The vorticity images convey information about the formation and dynamics of the vortices in the flow, which are intrinsically related to shear stress. We recall that Cheng et al. [25] showed the importance of low and oscillatory blood flow shear stress on the formation of atherosclerotic plaque; however, to the best of our knowledge, no work has been done to correlate the influence of shear stress on the endocardium and the cardiac valves. While the levelset approach as presented here would encounter difficulties in resolving the small geometric detail like the trabeculae carnae or cordae tendinae, if combined with the immersed boundary method it may be usable for evaluating valvular shear stresses.
5.1. Simulated cardiac cycle events
The simulation domain was discretized using a 128^{3} regular grid, which corresponds to a physical resolution of 1 mm, while the time step was 0.001 s. To ensure stability, we used subcycling whenever the maximum velocity would be greater than dx/0.001, i.e. the Courant–Friedrichs–Lewy (CFL) number max(u)*0.001/dx became greater than 1. The heart bounding box occupied about 95 per cent of the domain.
5.1.1. Systole
The particular heart we are focusing on has a systole of 290 ms in a cardiac cycle of 925 ms. The initial flow conditions inherited from the diastole and the isovolumic contraction (IVC) phases have the remnants of the toroidal vortex from the late mitral filling as the main feature (figure 5a). This feature dominates the lower ventricular pattern in the early systole (ES), such that the apical blood pool is mostly recirculated by the vortex. Interestingly, the initial conditions also feature a counterclockwise rotation, as seen from the apex, which reverses into a clockwise one during the systole.
The midsystole (figure 5b) features a strong aortic flux, often with vortical strands aligned with the aortic axis, which guide the righthanded helix motion of the blood into the aorta. This rotation of the blood continues downstream in the lower ascending aorta, and is a wellknown feature of the healthy aortic flow.
Our computational results are based on meshes obtained from CT data, hence not including torsional motion. They support the conclusion that the aortic blood righthanded helical rotation is mainly determined by the geometric disposition of the aortic longitudinal axis with respect to the aortic valve base plane and the aortic valve geometry itself, rather than the LV torsional motion. In endsystole (figure 5c), the vortical strand is weaker, and the fluid particles entering the aorta describe wider spiralling paths.
The flux is strong as well in the right heart's early and midsystole, with similar features as above, however, a regurgitant flux (measured and displayed in figure 8b) is visible at the end of the systolic phase, owing to a regurgitating pulmonary valve.
As can be seen in figure 6, our computations recover the known formation of the flow patterns distal to the aortic valve, namely the vortices forming behind each leaflet in the sinus region.
5.1.2. Diastole
The early diastolic (ED) flux begins with the formation of a asymmetric toroidal vortex (figure 7a) as the mitral valve opens, that further travels towards the apex at an angle of about 30° with respect to the axis linking the apex and the centre of the mitral ring (AMR). Upon hitting the posterior wall of the LV, the vortex suffers some dissipation and a change to its axis of about 75° with respect to the AMR. This evolves into a vortical pattern (figure 7b) that sweeps the apex with a large vortex with horizontal axis of rotation, which ends with a small vortex rotating in the opposite direction, located just below the aortic valve. The formation of this small vortex is enhanced by the anterior leaflet of the (open) mitral valve. This vortical pattern weakens gradually during middiastole, becoming barely visible later on. The late mitral filling (figure 7c) is successfully captured. It creates an extra toroidal vortex that travels downwards as described above (ED) and happens concurrently with a small pulmonary vein reflux, which is a normal event.
The flux is qualitatively similar in the right heart. We captured the reduced tricuspid flow curve diastasis compared with the mitral flow curve, as shown in figure 8b.
Our computations support the conclusion that tip vortices are being shed from the mitral valve during diastole. This is in partial disagreement with the view that the leaflets align themselves to the flow and therefore shed no vortices and do not steer the flow [4]. This difference is possibly owing to our enhanced heart model that includes accurate mitral leaflets. Our flow indicates that the mitral valve plays both a steering and a vortex shedding role. See, for example, figure 12b, where the mitral toroidal vortex connects to a vortex sheet that covers most of the upper surface of the mitral valve. In this simulation, the vortex sheet appears immediately after the start of the diastole, and separates at the mitral edge to create the toroidal vortex.
5.2. Quantitative analysis and discussion
For a quantitative evaluation, we estimated the blood flux across the valves. Typical phasecontrast magnetic resonance imaging (PC MRI) flow quantification protocols acquire blood flow across time using MRI planes, aligned with the anatomy by the operator (e.g. [26]). In order to perform a comparison with literature reported observations, we simulated this protocol and quantified the timedependentcomputed blood flow in the same fashion using a plane aligned with the valves as done for PC–MRI. The integral of normal velocity across this plane was computed and the obtained curves for the examined case are depicted in figure 8.
The flow is qualitatively similar to the standard flow curves for a normal heart [26], with the exception of the pulmonary arterial flow, which displays the regurgitation we mentioned before. Figures 9 and 10 show the two locations of the slice positions that were used to measure the flow. Note that the clinical practice for measuring average or peak velocities (for example, in ultrasound) uses slightly different measurement sites (the aortic site being slightly before the valve region, and the mitral one slightly below the mitral valve). We are, however, measuring blood flux (the integral of normal velocity across a given surface), therefore, we chose locations that are meaningfully displaced from (but still close to) clinical measurement sites, such that they

— separate the heart model into two topologically disconnected regions along the direction of the flow;

— do not cross the valve regions, which would possibly pollute (or at least make more complex) the flux computations.
Since the fluxes are measured (just) outside the LV, they may include volumetric change components owing to the radial motion of the aorta or of the atrium. While rather small, these components do show in figure 8 as regions of concurrent aortic–mitral (or pulmonarytricuspid) flow. For the same reason, the IP are not exactly quantifiable with this protocol, and therefore are less visible in figure 8 (especially the isovolumic relaxation in the LV).
The velocity patterns across the measurement planes are very complex as can be seen, for example, in figures 9 and 10, and this emphasizes the importance of using correct geometric information, especially for the valves, so that one can minimize the use of incomplete models for inflow or outflow. On the other hand, the flow curves from figure 8 would be unaffected by valve shape variation, with the exception of the valve opening and closing moments. While our paper shows the importance of the valves on the velocity and vorticity haemodynamic patterns in the heart (see also §5.3 on pathological cases), a more thorough study on the sensitivity of haemodynamics on the valve shape variation is needed.
As our study uses fourdimensional CT data, the contraction pattern of the ventricle is missing the torsional component, which can be as high as 15°. While we believe that such motion has little influence on the bulk features of the flow, we also note that—as emphasized in §4.1—the nonNewtonian nature of the blood may increase this influence. At the present time, there is no study of such influences, and a more thorough study is required. One way to include a torsional component is to consider mechanical models which include fibredirection.
5.3. Influence of pathologies—atrial modelling
Besides the pulmonary valve regurgitation, we evaluated other pathologies, such as a heart with a bicuspid regurgitant aortic valve and a pathological mitral valve, suffering from both stenosis and regurgitation. The blood flow pattern (figure 11) is very different than the normal pattern in the lower ascending aorta (figure 11h,i). The systolic jet, normally directed along the centreline of the aorta, is being deflected towards the aortic wall, increasing the local wall stress. This could explain why the patients with bileaflet aortic valve develop aortic root dilation. Furthermore, the pathological mitral valve, suffering from both stenosis and regurgitation, directs the flow towards the posterior LV wall (figure 11a,b,f). Regurgitations of the aortic and mitral valves can be observed during systole and diastole, respectively (figure 11g,j). This computation was initially reported in Mihalef et al. [9], and emphasizes the importance of the comprehensive cardiac model that we are proposing (especially the automatic and accurate detection of the valve motion) to the haemodynamic simulations.
We conclude this section with a few remarks based on our experience regarding the modelling of the atrium. Using the framework presented in §4, we performed simulations using models of the left heart with and without the atrium geometry (all the other parameters being identical). The systolic flow was almost identical, as expected (figure 12a,e). Schenkel et al. [4] and Long et al. [27] emphasized how big an influence the mitral boundary conditions had for their simulations. Our simulations showed that the main source of the influence is due to the vorticity generated mainly from the pulmonary veins, and transported downstream into the LV. This is clearly visible in figure 12b,f and c,g, which also shows that without the LAtoLV vorticity transport, the simulated flow inside the LV may not differ significantly between models with or without the LA. This is quite interesting, as it hints at a possible solution for modelling realistic boundary conditions without the LA, by including an appropriate model for the vorticity generation inside the LA, and its transport into the LV.
6. Conclusions
In this paper, we presented a comprehensive patientspecific modelling and computation of the human heart dynamics and haemodynamics. Geometric and kinematic parameters for the left and right ventricles, atria, cardiac veins, arteries and also cardiac valves, are estimated from fourdimensional cardiac CT scans using a robust learningbased algorithm. The estimated parameters are employed to enforce oneway momentum transfer to blood, using a levelsetbased Navier–Stokes solver. The final computations are demonstrated to provide flux curves that are qualitatively in agreement with standard literature PC–MRI measurements from comparable patients.
We are currently extending our simulation framework to include fibredirection and torsional motion, impedance boundary conditions at the outlet (the aorta) or inlets (pulmonary veins) and fully coupled FSI simulations by taking into account the twoway momentum transfer between the fluid domain and the solid wall.
Acknowledgements
This work has been partially funded by the European Commission through the project SimeChild (FP7ICT20094, no. 248421). We would like to thank Tommaso Mansi for help with automating certain visualization procedures and helpful feedback, and to the Interface Focus reviewers.
Footnotes
One contribution to a Theme Issue ‘The virtual physiological human’.
 Received November 16, 2010.
 Accepted February 25, 2011.
 This Journal is © 2011 The Royal Society