|
 
 |
| REVIEW ARTICLE |
|
| Year : 2019 | Volume
: 5
| Issue : 4 | Page : 141-153 |
|
Input data for computational models of blood flows
Marc Thiriet
Laboratory Jacques-Louis Lions, Sorbonne University, Pierre and Marie Curie Campus, F-75005, Paris, France
| Date of Submission | 12-Nov-2019 |
| Date of Decision | 10-Jan-2020 |
| Date of Acceptance | 10-Jan-2020 |
| Date of Web Publication | 13-Apr-2020 |
Correspondence Address:
Marc Thiriet
Laboratory Jacques-Louis Lions, Sorbonne University, F-75005, Paris
France
 Source of Support: None, Conflict of Interest: None  | Check |
DOI: 10.4103/digm.digm_24_19
Biomechanical models of blood flows in the last century were devoted to the handling of flow pattern in simple models more or less representative of rigid and deformable segments of the vasculature with bends and branchings, in normal and pathological conditions. With the development of medical three-dimensional (3D) imaging techniques, availability of 3D reconstruction software, improved fluid–structure interaction methods, and high-performance computing, computational models of blood flows are now aimed at guiding medical decision in the framework of precision medicine. However, most often, input data remain unavailable for a given patient, especially those related to blood and vascular wall rheology. This review summarizes the main blood flow features and yields the state-of-the-art.
Keywords: Blood flows, computational models, input data
How to cite this article:
Thiriet M. Input data for computational models of blood flows. Digit Med 2019;5:141-53 |
| Introduction | |  |
Blood is a flowing tissue, its matrix being plasma, in which are immersed blood cells (leukocytes) and pseudo-cells. Leukocytes are also called white blood cells, inflammatory cells, and immunocytes. Platelets are cell fragments containing storage vesicles, which participate in blood coagulation and inflammation. Red blood capsules (97% of “blood cells”) correspond mainly to a solution of hemoglobin (oxygen carrier) enclosed in a deformable plasma membrane.
The heart is constituted of the serial and apposed left and right pumps, each made up from an atrium and a ventricle to coordinate blood entry into and exit from these two pumps and to cope with upstream and downstream very low and more or less high pressures, respectively, for blood aspiration and ejection. In particular, blood is ejected simultaneously from the apposed left and right ventricles through the serial systemic and pulmonary vasculature. As the pressure is higher in the systemic circulation than in the pulmonary circulation, the myocardium of the left ventricle is more developed than the right ventricle one. Consequently, its contraction–relaxation cycle serves to describe the cardiac cycle.
Both circulatory circuits, systemic and pulmonary, are successively constituted from the ventricle of the corresponding pump to the atrium of the apposed pump of arteries, arterioles, capillaries, venules, and veins. Arterioles, capillaries, and venules form the microcirculation, whereas large, mid-size, and small arteries and veins are compartments of the macrocirculation.
Blood flows at the ventricular exit, that is, in the entry segment of the aorta and pulmonary trunk, according to a starting–stopping pattern. This pattern evolves to a pulsatile flow owing to proximal deformable elastic arteries, which store a fraction of blood ejection volume and restitute it during the following diastole using elastic recoil. On the other hand, flow resistance is set by small distal muscular arteries that exhibit a slight constriction state in the absence of stimuli. Vasodilation or constriction determines the local vascular lumen size and hence resistance. The vasomotor tone is regulated by the local stress field detected by cells of the vascular wall as well as autacoids, circulating messengers (e.g., hormones and remotely acting regulators such as angiotensin-2), and neurotransmitters released by nerve endings.
As any physiological system, the cardiovascular apparatus is characterized by a set of major properties: (1) diversity, that is, a huge between-subject variability in architecture (anatomy) of the vascular circuit, which explains the need of image-based three-dimensional (3D) reconstruction of the network for precision medicine (i.e., personalized or subject-specific investigation); (2) variability, that is, adaptation to environmental conditions, means that images acquired at a given time represented not only a model of the reality but also a frozen structure, which neglects influences acting on a living organ; (3) complicated structure and function are regulated both remotely by the fast-operating nervous as well as long-term–acting endocrine systems and locally by a locoregional control exerted by hemodynamic stresses, secreted autacoids, and metabolism of perfused organs; (4) complexity is illustrated by the heart pump that has a chaotic behavior in the deterministic framework, which enables it to rapidly responds to the environmental stimuli. A constant cardiac frequency yields a bad prognosis similar to strong anomalies in the generation and propagation of the cardiac action potential through the nodal tissue, which can engender an anarchic behavior of cardiomyocytes. This property that is beneficial with respect to the system functioning becomes a drawback in signal and image processing that involves averaging to enhance to signal-to-noise ratio.
As any physiological system, the cardiovascular apparatus is functioning with various length and time scales. Length scales are related to mechanisms that govern function of the cardiovascular system and its response to a changing environment, from: (1) cell signaling implicating messengers, receptors, and effectors (nm), such as molecules involved in the excitation–contraction coupling in cardiomycytes or in the regulation of the vasomotor tone by the couple formed by endotheliocytes and smooth myocytes; (2) adapting cells and their organelles (μm); (3) cell clusters organized in a tissue within an extracellular matrix (mm); and (4) organ, the cardiac pump or a vascular segment (cm). Multiscale modeling is aimed at coupling these different length scales, whereas the major objective of a multilevel modeling is to model the entire vascular circuit coupling 3D, 1D, and 0D model of the bloodstream.
While 3D models describe the field of hemodynamic variables (flow velocity and stress) in a set of nodes within a mesh obtained from discretization of a continuum, the distributed (1D) and lumped parameter models (i.e., electrical analogs) are strongly simplified versions of the reality. In a distributed parameter model, a vascular segment is assumed to be a succession of infinitesimally long slices, in which hemodynamic variables are only computed at the vessel axis, the axial tension between slices being neglected. In a lumped parameter model, the simplification increases a step further as hemodynamic variables are computed in a single node assumed to represent a given vessel between two branching points.
The time scales range from the fast adaptation relying on signaling axes associated with stored molecule release (ms to mn) up to gene transcription (h) and adaptive and adverse wall remodeling (e.g., cardiac response to regular exercise and hypertension; i.e., d, wk, and mo; [Figure 1]). For example, mechanosensitive (hemodynamic stress-gated) ion channel opens almost immediately (within ms) upon load bearing and elevation of cytosolic calcium (Ca2+) concentration and release of nitric oxide (NO) happens after a very short time (within s). In addition, the cardiovascular apparatus is subjected to the circadian cycle (i.e., a major day–night rhythm related to the sleep–wake cycle and accessory signals, such as food intake and body temperature). | Figure 1: Aorta branching and ligamentum arteriosum: a small side branch (upper right edge) and the closed entry of the large fetal ductus arteriosus (opposite edge). The wall thickness is not uniform, due to wall remodeling in this 3-month old calf by flow stresses, according to the nonplanar strong curvature of the aortic arch, using the locus of the ligamentum arteriosum as a landmark
Click here to view |
As hemodynamic factors are implicated in vascular diseases via the local stress field (pressure and wall shear stress), impingement force, and residence time of conveyed molecules, mechanical simulations complement the medical check-up.
Models usually rely on strong simplifications of the reality. The first step consists of collecting data (definition stage). Data relevant to the problem of interest are then selected (representation stage) as well as the theoretical framework and hence the governing equations. An appropriate numerical method is chosen while numerical tests carried out to display the fields of hemodynamic variables, to analyze results sensitivity to modifications of input parameters, to predict effects of a given set of conditions, that is, the output or a given set of input data, and to eventually propose changes in experimental conditions and optimized the medical procedure, surgical reconstruction, and/or design of implantable medical devices. In addition, inverse problems can be solved to assess physiological parameters that cannot be directly measured. Supposing that the numerical technique was verified, a validation stage is necessary due to consequences of generality loss associated with the model. Numerical results are thus compared to available observation data, either obtained by measurements in phantom or, better,in vivo using functional imaging techniques (velocimettry) coupled to anatomical imaging. Eventual iteration stages are aimed at refining the model.
Although the complete physiological system, the functioning of which depends on local and global factors, that is, intrinsic and environmental parameters, is more than the sum of its parts, reductionism remains a useful strategy of explore a complex reality, not only because the simple is solved before investigating the entire system, but also because a basic knowledge acquired on system parts is needed before developing a model of a complexity, the biological aspects of which remain more or less well handled.
Nevertheless, integrative models have been implemented and optimized. For example, a software platform couples the cardiac action potential (trigger), myocardial contraction, elevation of pressure to reach and overcome the outlet pressure load, and blood ejection,[1],[2],[3],[4] which is followed by relaxation and intraventricular pressure decay to enable ventricular filling. Integrative models incorporate behaviors at various length and time scales to efficiently describe the structure–function relation of physiological systems. Phenomenological models are yet employed to couple nano- and micro-scale events to the meso- and macro-scale behavior.
Whatever the focus, either on the cardiac or vascular wall, or the flowing blood, which can be possibly coupled in some studies, the theoretical framework is given by continuum mechanics. As blood flows most often in deformable vessels, except at least the cerebral superficial arteries and veins that run between the brain, an incompressible organ, and the skull, the fluid dynamics is strongly coupled with the solid mechanics. However, to conceive proper software is not trivial as the coupled media have similar physical properties. When dealing with walls of large arteries or veins, solid is usually assumed to be a thin-walled hyperelastic body that bears large deformation and obeys to the constitutive law given by an anisotropic Mooney–Rivlin material. The state-of-the-art in fluid–structure interaction remains at a preliminary stage, as the multilayed wall, each of its three tunicae being made up from a composite material, is modeled by a shell. In any case, even an adequate model of the wall of large vessels requires rheological data, which vary according not only to the spatial coordinate (i.e., longitudinal, azimuthal, and radial) but also to the wall layer. Their values in situ are not known. Any measurements carried out on isolated samples can provide values that differ from several orders of magnitude from those in vivo.
While arteries are distensible and dilate when crossed by the propagating pressure wave, veins are both distensible and collapsible. Veins constitute the blood storage compartment, as up to 70% of blood volume in the venous collector. However, veins may experience collapse, such as superficial veins of the inferior limbs, which are the site of dysfunctional valves. Stocking for varicose veins redirects venous blood flow from superficial to deep veins through valved perforating veins. Compression exerted of thin-walled superficial veins causes their collapse, but not that of deep veins embedded in tissues, which thus behave a thick-walled ducts. Unlike deformable pipes of neutral circular cross-section (when the transmural pressure equals zero), veins are elliptical. This geometrical feature avoids not only the occurrence of a weakly negative buckling pressure but also the bifurcation problem. During collapse, the number of lobes of a thin-walled tube of neutral circular cross-section depends on its geometry and rheology. Even for a tube with a weak ellipticity, the collapse is governed by a two-lobe collapsing mode.[5]
| Main Blood Flow Features | | |
In general, large vessels have a heterogeneous cross-section along their length between two branching stations, not only due to taper of long vessels such as the aorta but also due to prints of adjoining organs, even in the absence of compression. Their wall thickness in any cross-section is not uniform as wall adapts to the sensed local stress field [Figure 1].
Three-dimensional developing flow
The vasculature is characterized by a succession of bends of variable nonplanar curvature and branchings, either quasi-symmetrical (bifurcations) or, in general, asymmetrical. On the other hand, the venous network consists of junctions of draining venules and small, mid-size, and large veins. Nevertheless, the arterial bed possesses at least one junction, as the basilar trunk arises from the merging of the two vertebral arteries.
The vascular circuit is thus composed of various types of geometrical singularities: (1) curvatures in all directions that change continuously; (2) huge number of lateral branches, junctions, and terminal bi- and even tri-furcations; and (3) convergent and divergent segments, such as taper of a long vessel after successive branchings and prints on vessels imposed by apposed organs in normal conditions as well as in- and ex-trinsic pathological stenoses.
Vessel axis curvature causes helical fluid movement. Flow variables thus evolve in 3D fields confined by conduit walls. In addition, axes of some vessels such as coronary arteries and veins that travel in the epicardium over the heart surface bear large deformation, as they follow the motion of the cardiac wall during the cardiac cycle.
Any flow entering in a pipe experiences a resistance from the wall. The strongest effect of viscosity is exerted in the fluid boundary layer in the near-wall region of the vessel lumen. In the entry, or entrance, length (Le) of the vessel, the momentum diffuses toward the vessel axis and the boundary layer grows. The boundary layer thickness is the distance from the wall to a point where the flow velocity becomes equal to the core-flow velocity, in which the fluid is assumed to behave as an inviscid material.
A fully developed flow that appears when the boundary layer reaches the vessel axis in a straight pipe is characterized by an axially invariant velocity cross-distribution and a linear axial pressure drop (constant pressure gradient).[6] In any case, in a deformable bend, the short distance between two branch points prevents the appearance of a fully developed flow.
As blood flows in entrance segments, the vessel length being lower than the entry length for the local Reynolds number, the pressure decays nonlinearly.
Time-dependent flow
Blood circulation is characterized by a bulk unidirectional transport during most of the cardiac cycle in arteries. However, bidirectional flow, that is, blood flows in the streamwise direction in the core lumen and in the opposite direction in the near-wall region, the thickness of which is not uniform in a cross-section of a bend, appears in elastic arteries (close to the heart) during diastole, except in carotid arteries of the cerebral circulation. The carotid arteries do not have the same rheology as the homologous subclavian arteries. Moreover, the blood pulsatile flow reverses its direction throughout the entire cross-section toward the heart in some proximal arteries, such as the aorta. This flow behavior enables an efficient coronary perfusion during diastole, when the myocardium is relaxed, which is stronger than that during systole through the transiently contracted myocardium, the coronary ostia being partly hidden by the aortic valve leaflets.
Deformable wall
Physiological fluids such as blood are transported in deformable conduits. Distensibility of large proximal arteries (e.g., aorta and common carotid and subclavian arteries as well as large pulmonary arteries), which is linked to elastin-enriched walls, damps flow fluctuations during the cardiac cycle at the ventricle exit.
During the contraction of the ventricular myocardium, the ejected systolic blood bolus is not only partly transmitted in the perfusion network but also partly stored in large elastic arteries close to the heart exit. During the following ventricular myocardial relaxation, the stored blood is restituted and converts the starting–stopping flow downstream from the ventriculoarterial valves into a pulsatile stream in downstream arteries.
The locally and transiently stored blood flow during the systole is transferred in both up-and down-stream directions during the following diastole. Therefore, it maintains organ perfusion during diastole. In particular, upstream blood flux permits coronary irrigation, hence heart perfusion.
The Windkessel effect (from German Windkessel: air chamber) is associated with the capacitance of the large proximal elastic arteries. It is related to arterial dilation during systole, as blood pressure rises, and subsequent small constriction during diastole. The systolic storage of a stroke volume fraction is restituted by elastic recoil during the following diastole, transforming the starting–stopping flow at the heart exit into a pulsatile flow in arteries.
Any artery transiently dilates during the passage of the traveling pressure wave. Veins also dilate; they constitute the blood storage compartment. Veins of inferior limbs experience a partial collapse during contraction of the skeletal muscles between which they run, especially during walking.
Any tube of unstressed circular cross-section changes its transverse configuration when the transmural pressure (p; internal minus external pressure) is lower than a buckling pressure, which is slightly negative. Moreover, the number of lobes that appear after contact between opposite walls depends on the duct geometry and rheology. The lobes are the open part of the collapsed tube lumen that is associated with symmetry axes.
On the other hand, veins have a slightly elliptical unstressed cross-section (A). This geometrical feature avoids the bifurcation problem. Veins collapse according to a two-lobe mode. Furthermore, the cross-sectional area decreases abruptly as soon as the transmural pressure becomes negative.
Three major values of the transmural pressures characterize the collapse: (1) the ovalization pressure, for which the curvature at the points of the mid-line at the top and bottom faces is locally equal to zero; (2) the point-contact pressure, at which the opposite sides touch for the first time, the contact reaction inducing a discontinuity in the first derivative of the A (p) function; and (3) the line-contact pressure, when the curvature at the contact point equals zero, the contact reaction provoking a discontinuity in the second derivative of the A (p) function.
The collapse mode 1 corresponds to the collapse before contact, that is, when the transmural pressure is slightly negative and the cross-sectional area decreases sharply, the vessel being then characterized by a high compliance. Mode 2 is characterized by a contact at a single point, the curvature at this contact point decaying from a finite value down to zero. Mode 3 is defined by a contact on a line segment, which lengthens with the falling transmural pressure.
As communicating veins between the superficial and deep venous limb circuits connect vein edges, they remain irrigated by the drainage lobes of collapsed veins. However, vein collapse may not reach the contact configuration, that is, when opposite walls come into contact at one point with a local finite and then infinite curvature radius. Such an event is associated with a discontinuity of the transmural pressure–cross-sectional area relation as pointed out previously.
Flow pattern
Flow pattern is related to flow stability. In other words, flow pattern relies on the evolution of disturbances. Instabilities do not always engender a turbulent flow. When stabilizing viscous forces are strong, they damp flow perturbations. On the other hand, flow instabilities grow with time when inertia forces are very large.
The ratio of inertia forces to viscous forces corresponds to a dimensionless flow governing parameter, the so-called Reynolds number (Re). Hence, Re measures the strength of the nonlinear convection term with respect to the linear viscous term, the former amplifying flow instabilities. However, flow stability depends not only on Re but also on wavelength and propagation speed as well as evolution of flow perturbations.
In general, decay and growth of flow perturbations were studied when their amplitude was infinitesimal rather than finite. In addition, most investigations deal with infinitely long, uniform, straight, and rigid pipes that convey a steady flow. In the so-called laminar regimen, when flow is steady or not, when Re is smaller than a critical value (Recr), fluid moves through a straight duct in a somewhat orderly fashion, that is, in sheared laminae.
When Re increases Re > Recr, some regions of the flow become chaotic. Flow evolves toward a transitional regimen. Turbulence is intermittent and alternates with quiescent laminar regions. This transitional regimen may evolve to a fully turbulent state, a random motion being observed everywhere within the flow. Transition between laminar and turbulent pattern is initiated by a primary instability that produces secondary and generally 3D motions that become unstable.
In flow experiments, small perturbations that define the rate of residual turbulence are, in general, damped by fluid viscosity linked to momentum transfer between neighboring fluid elements. Whatever the flow experiment care, when the transition threshold is overcome, successive tiny perturbations are amplified and flow is disorganized.
In a long straight pipe of uniform circular cross section and smooth rigid wall, O. Reynolds observed a threshold above which inertia forces become predominant and flow becomes disorganized. This threshold was later named the critical Reynolds number. The usual value found in any textbook ranges from 2000 to 2300. However, it is impossible to accurately determine a value above which the laminar pattern disappears, because, in carefully performed experiments in straight pipes, Recr can reach an order of magnitude of 105. Moreover, values for physiological flows, that is, for oscillatory and pulsatile flows in deformable curved vessels, remain unknown.
Laminar flow corresponds to fluid particle motions in adjacent layers, usually when inertia forces remain less than three order of magnitude greater than viscous forces. When viscous effects are strong enough, the internal fluid coherence is preserved and fluid moves by layers. Poiseuille flow in a straight duct of circular cross-section, which is the first analytical solution of the Navier–Stokes equations (describing the fluid flow in a given domain with given initial and boundary conditions), is displayed by fluid particle motion in concentric laminae. Couette flow in a gap between two concentric cylinders, a second analytical solution of the Navier–Stokes equations, is depicted by fluid particle movement in parallel laminae.
The knowledge of a fluid particle motion, or a buoyant-seeded solid particle used for visualization at a given flow point at a given instant, enables prediction of the movement at any subsequent time, as the particle trajectories are similar. On the other hand, when Re > Recr, the correlation between particle trajectories crossing the same point at different instants is lost, especially when the time interval between two successive observations increases. The motion becomes unpredictable. A tiny gap in initial conditions provokes catastrophic consequences in long-term flow evolution.
The term laminar flow is a misnomer as it can be associated with predictable flow separation, occurrence of eddies, and vortex shedding due to geometrical singularities and/or immersed bodies, without loss in correlation between fluid particle motion. Errors on initial and terminal assessments have the same order of magnitude, without growth. In other words, imprecision on initial conditions does not prevent prediction.
Fluid swirling happens when pressure forces on the opposite fluid particle faces cannot overcome the torque engendered by the shear, that is, the vorticity of the fluid particle is unbalanced.
The value of the Reynolds number should be appropriate to flow conditions (e.g., time-dependent flow), vessel geometry (e.g., curvature), and wall rheology (rigid vs. deformable). The Reynolds number is the ratio of a product of the flow length and velocity scale to the fluid kinematic viscosity. When dealing with flow perturbations, the relevant length scale is the boundary layer thickness and not the vessel radius (i.e., the largest size of the boundary layer thickness of a fully developed flow in a straight duct).
In a stable flow, any small disturbance is damped. Strong local pressure gradients, velocity profile inflections, and shear reversals affect the flow stability. Flow instabilities can be related to vortex production; vortices arise from high-vorticity regions. Sheared vortex can then be stretched and rotate. Self-modulated perturbations can lead to evolving small-scale structures. In an intermittently disturbed flow, transient instabilities decay as they propagate downstream, due to the dissipative action of viscous forces. Amplification of vorticity perturbations in flow can lead to a transitional flow.
Periodic flows are transiently stable if a disturbance grows during a part of the cycle and decays. On the other hand, periodic flows are unstable when the disturbance grows during each period. Periodic flow patterns can be classified into four main types: (1) laminar; (2) disturbed with small amplitude perturbations; (3) intermittently turbulent flow in which high-frequency velocity fluctuations occur at the beginning of the deceleration phase, increase, and dissipate before or during the subsequent acceleration phase; and (4) fully developed turbulent flow, high-frequency velocity fluctuations existing during the whole flow cycle.
As flow disturbance can only grow during the decelerating phase, pulsatile flows are more stable than steady flows. In addition, many flow features affect the instability growth, such as vessel axis curvature for a given perturbation amplitude, wall deformability, and harmonic content of the pressure wave.
Living, sensing, and reacting wall
Large arteries and veins possess a multilayered wall that is composed of an inner intima, a media, and an outer irrigated and innervated adventitia. Each wall layer is made of a composite material. The adventitia contains a perfusing microvasculature (vasa vasorum). The intima is made of a single, wetted layer of endotheliocytes over a thin sheet of connective tissue. The media consists of sublayers of smooth myocytes and elastic fibers, the ratio of which determine the nature of the vessel. Elastic and muscular arteries are enriched of elastic fibers and smooth myocytes, respectively.
Contraction or relaxation of smooth myocytes determines the luminal caliber. In normal conditions, arteries are slightly constricted. This state, especially in small arteries and arterioles, determines peripheral vascular resistance.
Mural cells, endotheliocytes and smooth myocytes, sense the local stress field and react. Blood pressure (stress component normal to the vessel wall) engenders azimuthal and axial tension. Endotheliocytes are also subjected to the friction of flowing blood at the wall–blood interface (tangential stress component), the wall shear stress. Although endothelial and smooth muscle cells detect a priori both azimuthal and longitudinal tension induced by the blood pressure, which fluctuates with a relatively large amplitude during the cardiac cycle, endotheliocytes react primarily to the local shear, which evolves with a relatively low amplitude.
The parietal stress field operates via mechanogated and sensitive proteins embedded with the plasma membrane of vascular cells, the former reacting directly and the former indirectly via proper mediators.
Smooth myocytes and endotheliocytes contribute to the regulation of the vascular caliber directly and via release of vasoactive substances (vasodilators such as nitric oxide and vasoconstrictors such as endothelin), respectively. Vasoactive substances control the contraction–relaxation state of medial smooth myocytes, participating in the local regulation of the vasomotor tone. The latter is also controlled by neurotransmitters secreted by endings of the nervous system as well as circulating hormones and other chemical messengers, such as products of the metabolism of perfused organs. The flow can then adapt to maximize the organ perfusion and minimize the heart afterload.
Blood rheology
RBCs represent 97% of blood cells and thus play a major role on the rheological blood behavior. Blood can be indeed considered as a concentrated RBC suspension in a Newtonian solution (plasma), which is composed of ions and macromolecules interacting between them and bridging RBCs.[7]
In the macrovasculature, the ratio between the vessel bore and cell size is such that blood is considered as a continuous homogeneous medium. In the microvasculature, the flow length scale is small and blood is a heterogeneous medium flowing in corrugated vessels, due to the presence of endothelial protrusions associated with cell nuclei. In capillaries, blood flows as a train of deformed RBCs that entrap boli of Newtonian plasma and are surrounded of a lubrication layer of marginal plasma.
The rheological properties of non-Newtonian blood are dictated by the flow-dependent evolution of the blood internal microstructure, i.e., the state of the flowing RBCs (possible aggregation and deformation). The interactions between the conveyed plasma molecules and RBCs indirectly govern blood rheological behavior.
Shear-step experiments (steady state after a short transient regimen) show that the blood has a shear-thinning behavior.[8] The shear-thinning behavior is exhibited by a pseudosigmoid relationship between the logarithm of the shear rate and the logarithm of the blood apparent viscosity. This behavior is explained by RBC aggregation at low-shear rates (high apparent viscosity values), whereas RBC deformation and orientation lead to a low viscosity plateau μ∞ at high-shear rates  (Newtonian behavior range). A huge between-subject variability exists in the  relationship. This variability points out the need for standardized procedure (temperature, pH, Ht, protein concentration, etc.). The  relationship depends on Ht.[9] The very-low-shear-rate viscosity (μ0; maximum of blood apparent viscosity), which is difficult to measure with available experimental devices, is about ten times greater than the high-shear-rate viscosity (μ∞). When both RBC aggregation and deformation are inhibited, the blood behaves almost like a Newtonian fluid.[10]
However, blood is a time-dependent shear-thinning fluid (thixotropy). Its behavior is affected by time-dependent, reversible changes in the blood structure. The thixotropic behavior is indeed explained by aggregation and disaggregation of RBC rouleaux. (1) At rest, blood is composed of a rouleaux network, that is, a medium that does not represent the flowing blood, in which sheared RBCs are isolated and oriented in the streamwise direction. (2) Low-shear rates deform and partially break the rouleaux network into smaller rouleaux networks and isolated rouleaux. (3) At mid-shear rates, only isolated rouleaux of various sizes is observed. (4) At high-shear rates, RBCs are scattered and oriented in the shear direction. A shear rate step induces a transient and stable strain regimen, with a viscoelastic, elastothixotropic, and Newtonian response at low  , mid  , and high shear level  , respectively.[11]
Rheological data, which are provided by experimental results obtained in steady-state conditions (blood behavior depends on the structure of the RBC population), are far from the physiological ones in the arteries.[12] Shear-step experiments do not match the blood flow dynamics. The flowing blood is characterized by a smaller convection time scale than the characteristic time of erythrocyte aggregation. In addition, erythrocytes expelled from the left ventricle, where they have been shaken during the isovolumetric contraction, do not have time to aggregate in the absence of stagnant blood regions. On the other hand, non-Newtonian fluid flow must be simulated in stagnant blood regions (aneurysmal cavities and flow separation downstream from severe stenoses), when the residence time of erythrocytes is greater than the aggregation time constant.
In the framework of continuum mechanics, the time-dependent shear-thinning behavior and loading history of blood are incorporated via an appropriate formulation of the constitutive law. Generalized Newtonian models used to mimic the shear-thinning blood behavior do not take into account the stress and strain history of blood and aggregation–disaggregation cycles of the red blood cells. They are not relevant in pulsatile flows.
In stagnant flow region, in particular, when an aneurysm develops in the arterial bed with a configuration that promotes a local flow stagnation, red blood cells can have time to aggregate and form rouleaux. These rouleauxs can be modeled as purely elastic dumbbells, the length of which changes as rouleaux aggregates and fragments. The resulting suspension of isolated erythrocytes and erythrocytic rouleaux oriented in the streamwise direction constitutes an incompressible, shear-thinning, viscoelastic, and thixotropic fluid. The set of equations, which is similar to an Oldroyd model, to solve becomes:[13]
Where E represents the contribution of the elastic dumbbells to the total Cauchy stress,  represents the average size of rouleaux, μ = Ht (BoT + k)τ represents the polymer viscosity (Ht: Red blood cell concentration, Bo: Boltzmann constant, T: Absolute temperature,  : Relaxation time that depends on erythrocyte aggregation and fragmentation rates and relaxation time of a single particle),  represents the erythrocytic rouleaux fragmentation rate, and  represents the value of average aggregate size N in a simple shear, steady flow with shear rate  , where A (D) is the aggregation rate and τ∞the relaxation time when the shear rate tends to infinity. The ratio  .
In the elastic stress component due to aggregates of all sizes (various aggregated cell number [N]), the relaxation time  for aggregates of size Ni is replaced by an average value  .
| Physical Modeling of a Flowing Fluid | | |
Fluid flows are governed by mass, momentum, and energy balance principles, which are expressed by partial differential equations (PDE). The main fluid variables, the velocity vector v(x, t) and stress tensor C(x, t), use the Eulerian formulation (x: Eulerian position; t: Time). The Navier–Stokes equations that govern Newtonian fluid flow are nonlinear PDEs that result from the balance of forces in isothermal flows of an incompressible fluid in pipes.
The mathematical formulation of the equations governing movement of any infinitesimal control volume, the so-called fluid particle of mass density ρ, dynamic viscosity μ, and kinematic viscosity ν = μ/ρ, derives from loading analysis and mass transfer.
The first basic postulate states that matter is neither destroyed nor created, meaning that any finite volume of matter neither disappears nor becomes infinite. The second basic postulate states that matter is impenetrable; any matter element does not share the same location with another element (continuity axiom).
The following differential operators are used: (1) the gradient operator  (2) the divergence operator ∇; and (3) the Laplace operator  . The gradient ∇ p of scalarpis a vector of component  . The gradient ∇v of vector v is a second order tensor with component  .
Mass and momentum conservation are expressed by the following formulas:
where  (I: Metric tensor; D: Rate of deformation  for a Newtonian fluid). In isothermal conditions, the equation set (2) leads to the simplified form of the Navier–Stokes equation in the absence of body forces b:
When the inertia term (v⋅∇) v can be neglected, i.e., the diffusive term  , is predominant, such as in small arterioles, the momentum conservation equation is called the Stokes equation:
Main flow governing parameters
Acting forces encompass: (1) temporal inertia forces  V: flow velocity scale, ω: flow pulsation representing the flow time scale); (2) convective inertia forces  L: flow length scale); (3) viscous forces  (4) pressure forces; and (5) possible body forces, such as gravity  and magnetic field during nuclear magnetic resonance (MR) imaging.
Hence, the main flow governing parameters comprise: (1) Reynolds number (Re), that is, the CIF/VF ratio (UL/ν); (2) Stokes number (Sto), that is, the TIF/VF ratio (L(ω/ν)1/2); and (3) Strouhal number (St), that is, the TIF/CIF ratio (ωL/U). Their range of values within the macrocirculation is given in [Table 1]. The Stokes number is also commonly called Womersley number and sometimes the Witzig-Womersley number in the literature. The Stokes number also defines the time ratio  as well as the length ratio R/δS, where δS= (v/ω)1/2 is the Stokes boundary layer thickness. While the Rayleigh boundary layer deals with a flow over a flat plate which suddenly moves in its own plane, with a constant speed (i.e., a transient regimen), the Stokes boundary layer is linked to a harmonic motion of a flat plate in its own plane, with an angular frequency ω (i.e., periodic flow), which is more relevant to physiological flows, hence the preferred name “Stokes number.” | Table 1: Decreasing values of the main dimensionless flow governing parameters from the entry segment of the aorta to a small artery
Click here to view |
The Navier–Stokes equations suitable for a bend of uniform curvature in a single curvature plane using toroidal rather than Cartesian coordinates can be found in the literature.[14]
Additional dimensionless parameters can be deduced from similar analyses. The Dean number (De = [R/Rc]1/2 Re; Rc: curvature radius of the tube axis) for a laminar flow in curved vessels is the product of the square root of the vessel curvature ratio by the Reynolds number. It is usually calculated in simple bends of constant curvature, such as those used in experiments or simulation benchmarks. In image-based flow models, the vessel axis varies continually in every direction; consequently, the Dean number is not computed.
The modulation rate, or amplitude ratio, is used when the time-dependent component of the fluid flow is a sine wave of amplitude V∼ superimposed on a steady one  to assess the magnitude of the near-wall back flow during the bidirectional phase of the pulsatile flow. The lower the modulation rate, the smaller the retrograde flow. Blood flows are engendered by pluriharmonic pressure waveforms.
Dimensionless form of Navier–Stokes equations
The Navier–Stokes equations can be transformed into a dimensionless form using appropriate scales for the length (L; either the boundary layer thickness or duct hydraulic radius for fully-developed flow), time (T; either the flow cycle period [T] or inverse of the flow pulsation  , velocity (V; either the mean or maximal cross-sectional velocity), and the pressure (P).
The classical formulation of the dimensionless Navier–Stokes equations is hence obtained:
According to the predominant force type (convective inertial forces or viscous forces), the dimensionless coefficients related to the terms T1– T4 of the equation 5 are displayed in [Table 2]. In unsteady flows, the local inertia effects can be preponderant. | Table 2: Dimensionless coefficients  of the Navier-Stokes equation in the absence of body forces (equation5). Cp, Na, Re, Sto, and St are the pressure coefficient and Navier, Reynolds, Stokes, and Strouhal numbers  : flow transition)
Click here to view |
Both the Stokes and Strouhal number are important to assess the flow unsteadiness [Table 3]. In steady flows, the pressure gradient balances the friction (and gravity when it has significant effect, as in blood vessels in standing position). In quasi-steady flows, the rate of change in boundary conditions is so slow that the momentum has time to diffuse during the flow period. | Table 3: Dimensionless governing parameters of arterial blood flow (cardiac frequency 1 Hz, blood density of 1055 kg/m3, dynamic viscosity 3.5×10-3 Pa·s)
Click here to view |
Physiological boundary conditions
Boundary conditions, as well as initial conditions, are required to solve the set of PDEs. Numerical simulations of the flow in a segment of physiological bioconduit network require the specification of boundary data on artificial boundaries of the computational domain, which limit the explored vascular district. The boundary of the fluid domain (Ω) is usually partitioned into three main types of surfaces: vessel entry (Γ 1) and exit (Γ 2) cross-sections and wall (Γ 3). The classical no-slip condition is applied to the rigid vessel wall. On the other hand, in the case of a deformable wall, the moving boundary Γ3 is the fluid–structure interface. Inlet and outlet boundary conditions influence the fluid dynamics and must then be set away from the region of interest. Moreover, upstream and downstream effects caused by any geometrical singularities such as bends must be incorporated; hence, the 3D computational domain must contain proper boundary regions at some distance from the explored region to avoid numerical bias.
In the cardiovascular system, flow and pressure waves emanate from the heart and travel through the major arteries where they are damped, dispersed, and reflected due to changes in vessel caliber and wall properties as well as branchings. As a consequence, solutions to the governing equations of blood flow in large vessels depend on inlet and outlet impedances.
A time-dependent injection velocity vΓ(t) can be prescribed, at least, at the inlet. Most often, the spatial resolution ofin vivo velocity measurements is indeed not high enough to provide velocity distribution at vessel ends. Velocity conditions can be obtained from Fourier transforms ofin vivo ultrasound or MR flow signals. These measurements are commonly carried out more or less far from the computational domain entry and not necessarily in the same subject.
In the absence of measurements, the inlet velocity profile is given by the Womersley solution.[15] The Womersley solution corresponds to a Poiseuille flow superimposed to a harmonic flow, hence to a pulsatile fully-developed flow in a long, smooth, straight pipe of circular cross-section, a condition never encountered in the human body. On the other hand, time-dependent, uniform injection velocity is associated with an artificial high wall shear, hence vorticity, at the entrance station, that is susceptible to induce larger flow separation, if an adverse pressure gradient is set up.
Bulk interactions between the region of interest and network parts upstream and downstream from it are currently neglected. In fact, at the outlet cross-sections, the flow distribution and pressure field are not measured. The most common outflow boundary conditions for 3D flow simulations are prescribed traction and velocity profiles (The outlet sections must be perpendicular to the local vessel axis and be short straight pipe exits to avoid pressure cross-gradient):
where n and t are the local unit normal and tangent vectors and the outlet pressure is approximately equal to −f0. In the case of multiple duct exits, f0 influences the flow distribution among the model branches. These standard boundary conditions at intlet(s) and outlet(s) of a small explored region of the vasculature do not take into account the input and output impedances and then cannot predict accurately the flow behavior.
Artery walls can bear large displacements, as the arterial lumen can vary up to about 10% between the diastolic and systolic configuration. A multiphysics approach is then necessary to provide suitable boundary conditions both at the moving fluid–solid interface and domain inlets and outlets.
When physiological measurements cannot produce suitable boundary data, a mathematical description of the action of the reminder of the fluid circuit on the studied region can be carried out.[16] A multilevel modeling couples 3D flow models (detailed models relying on PDEs) to 0D (lumped-parameter models, or electrical analogs, described by ordinary differential equations). In the case of deformable bioconduits, 1D models are incorporated for adequate wave propagation. However, although it incorporates the vascular circuit, the multilevel modeling uses vessel parameters that are not related to the subject from which the 3D region of interest is explored.
| Blood Flow Modeling | | |
Numerical experiments of blood flows rely on a multimodeling process. (1) Medical imaging, from which originates the computational domain, is an anatomy modeling. (2) Facetization associated with the 3D reconstruction is a model of the actual geometry of explored vessels and architecture of the vascular circuit. (3) Spatial discretization yields a model of a continuum. (4) Simulations are based on physical, mathematical, and numerical approximations. Therefore, the simulation goal is to give qualitative information rather than quantitative one.
Modeling drawbacks comprise: (1) a fixed geometry and hemodynamical conditions, that is, a frozen normal or pathological state at rest in the lying position; (2) results that are more or less dependent on the image processing methodology mesh quality; (3) weak handling of the flow rate distribution in the explored circuit that can vary during distinct instants of the cardiac cycle, when it was measured; (4) boundary conditions, functional data are, most often, not obtained in the subject, the anatomy of whom is investigated, and either stress-free conditions are utilized at the exits, or simplified models are used, but rheological data are not those of the subject.
Three-dimensional reconstruction
Computerized medical imaging provides subject-specific 3D geometry of any body organ, in particular large vessels. The common two-step technique to create a mesh from imaging data consists first in segmenting the selected vessels and then using the segmentation surface to create the facetization (surface discretization). Most often, automatic mesh generators are used; they must be able to cope with surface that is frequently full of gaps, overlaps, and other imperfections. Surface meshes must not only be accurately approximate the vessel shape but also be suitable to computations, especially at vessel ends.
Mesh
Mesh element density should not influence numerical results. In addition, an optimal 3D mesh of the vascular circuit is refined near the vascular walls and coarsens in the core region around the vessel axis.
The surface discretization obtained from the 3D reconstruction needs further treatment to be suitable for numerical tests. Boundary conditions must be set at domain boundaries far from the exploration volume to incorporate both upstream and downstream effects of geometrical singularities on flow. Short straight ducts in the direction of the local axis terminating with a cross-section to avoid transverse pressure gradient are connected to every vessel end.
The methodology of surface reconstruction also affects numerical results. This effect is illustrated by arterial contours extracted from a set of saved scan slices with a between-slice distance of 3 mm focused on the common iliac artery with a side aneurysm and its bifurcation into two arteries [Figure 2]. The 3D surface is reconstructed using three different methods, a spline-fitting method associated with a CAD software (model M1) and two home-made software dedicated to surface reconstruction from polygonal contours.[17] The surface of model M2 is reconstructed from planar parallel contours using the “NUAGE” software, which is mainly based on Delaunay triangulations and Voronoi diagrams of points located in two parallel planes. The same set of input points is used for model M3. After a cubic spline fit of each contour associated with a smoothing, a new node set is defined by equally spaced points along the vessel contour. Then, two successive slices are projected orthogonally to the local axis in a same plane and a 2D-constrained Delaunay triangulation is built. The surface triangles are finally extracted by elevation of the two planes, the set of intrinsic axes of arteries giving the projection direction for the slice pair treatment, whereas in M2, a general projection direction normal to the slice planes is used. | Figure 2: Various reconstructuion techniques of a side aneurysm on the common iliac artery (from left to right, model M1, M2, and M3)
Click here to view |
The meshing techniques have their respective advantages and drawbacks. The model M1 possesses the best aneurysmal nucha edge that exhibits a slight vertex, which is observed in the control angiography (flat nucha in both M2 and M3). It has also a better aneurism height-to-width ratio. The model M2 displays a better ratio between the neck and sac width as well as a better aortic bifurcation apex. The model M3 reproduces with a better quality the nose front edge as well as arterial branches. However, the configuration of the computational domain affects the field of the hemodynamic quantities. The relative difference in pressure maximum is about 25% smaller in M3 than in M2. Nevertheless, this difference may not be significant in the context of sequential multimodeling involved in image acquisition and processing, meshing, and numerical approximations.
| Experimental Validation | | |
Validation of numerical results is illustrated by comparison with measurements of a steady laminar flow of incompressible Newtonian fluid in a model of the human common carotid artery and its main branches, especially the internal carotid artery (ICA) and external carotid artery (ECA). Steady flow in the carotid artery mimics the situation after implantation of a small continuous-flow ventricular assist device.
A P1–P1 stabilized SUPG/PSPG/CONT finite element method is used to solve the Navier–Stokes equations in the laminar pattern. Blood flows are simulated at local Reynolds number 320 (based on vessel radius and time-mean cross-sectional average velocity in pulsatile condition) for different but plausible outflow conditions, i.e., flow distribution between ICA and ECA. Simulation errors using a Zienkiewicz-Zhu, a posteriori error estimator, are computed as well as local gradients in dynamical quantities. Nine grids are used to ascertain mesh independence. A measured entry velocity distribution is imposed at the inlet of the computational domain. Either traction-free conditions or given values associated with lumped parameter models of the downstream vasculature are set at outlets.
A silicone rubber phantom of the computational domain has been built by rapid prototyping. Neutrally buoyant particles are seeded in a fluid, the refractive index of which matches that of the phantom wall. A time-resolved stereoscopic particle image velocimetry serves to measure the velocity field. Simulations agree with measurements [Figure 3].[18] | Figure 3: Numerical (left) and experimental velocity contours (right) in selected ross sections of an image-based carotid arterial circuit conveying a steady flow mimicking the presence of a ventricular assist pump
Click here to view |
The flow ratio is most often assumed to be equal to 7:3 (between internal and external branches) under normal conditions. However, with ventricular assist device (VAD), the flow ratio has not been assessed. Therefore, in the present work, the flow ratio ranges from 7:3 to 3:7. The outflow boundary conditions do not markedly affect the flow in the trunk upstream from the transition zone. With a flow ratio 7:3, the flow field at the bifurcation divides between the two branches without strong recirculation in the carotid sinus. With decreasing flow ratio, high-velocity isocontours in the transition zone moves toward the stem axis. With a flow ratio 3:7, the flow is deported into the external carotid increasing the strength of the separation region in the bifurcation segment. The maximum pressure point is located near to the maximum wall shear stress on the carotid bifurcation. These values increase in intensity as the flow ratio rises from 7:3 to 3:7. The arterial wall undergoes subsequent remodeling under continuous flow after VAD implantation.
| Concluding Remarks | | |
Blood flow in large vessels depends on time. Curvature and branchings are responsible for a developing 3D flow. Blood flows through deformable vessels, the rheology of which remains to be assessed in vivo. In the absence of instability growth, blood flow is laminar. In the absence of blood stagnant regions, blood behaves as a Newtonian fluid.
Boundary conditions at vessel inlets and outlets as well as flow distribution among branches can rely on functional imaging. When the vascular cross-section is too small with respect to the spatial resolution of the imaging apparatus, the vascular circuit can be modeled by a set of lumped parameter models. However, the associated parameters, which depends on the subject, especially in pathological conditions (e.g., serial arterial lesions such as stenoses), remain to be determined.
On the other hand, blood flows as a series of deformable cells between plasma boli surrounded by a lubrication plasma layer, at least in capillaries of caliber greater than 5μm. This particle flow, which ensures exchanges of gas, nutrients, and wastes, can be assumed to be steady. In capillaries, endothelial corrugations linked to the presence of cellular nuclei should not be neglected and wall law used for boundary conditions.
Because of the multimodeling methodology, computational models are aimed at providing qualitative data rather than quantitative ones. Nevertheless, these qualitative data can improve diagnosis, treatment planning, and prognosis, thereby assisting medical decision.
Financial support and sponsorship
Nil.
Conflicts of interest
There are no conflicts of interest.
| References | | |
| 1. | Bestel J, Clément F, Sorine M. A biomechanical model of muscle contraction. In: Niessen WJ, Viergever MA, editors. Medical Image Computing and Computer-Assisted Intervention (MICCAI'01), Lecture Notes in Computer Science (lNCS). Vol. 2208. Switzerland: Springer; 2001. p. 1159-61.  |
| 2. | Sainte-Marie J, Chapelle D, Cimrman R, Sorine M. Modeling and estimation of the cardiac electromechanical activity. Comput Struct 2006;84:1743-59. |
| 3. | dos Santos ND, Gerbeau JF, Bourgat JF. A partitioned fluid-structure algorithm for elastic thin valves with contact. Comput Methods Appl Mech Eng 2008;197:1750-61. |
| 4. | Smith N, de Vecchi A, McCormick M, Nordsletten D, Camara O, Frangi AF, et al. EuHeart: Personalized and integrated cardiac care using patient-specific cardiovascular modelling. Interface Focus 2011;1:349-64. |
| 5. | Thiriet M, Naili S, Langlet A, Ribreau C. Flow in thin-walled collapsible tubes. In: Leondes C, editors. Biomechanical Systems. Techniques and Applications: Biofluid Methods in Vascular and Pulmonary Systems. Boca Raton, Florida: CRC Press; 2001. |
| 6. | Thiriet M, Naili S, Ribreau C. Entry length and wall shear stress in uniformly collapsed veins. Comput Model Eng Sci 2003;4:473-88. |
| 7. | Thiriet M, Sheu TW, Garon A. Biofluid flow and heat transfer. In: Johnson RW, editor. Handbook of Fluid Dynamics. 2 nd ed. Boca Raton: CRC Press; 2016. |
| 8. | Chien S. Shear dependence of effective cell volume as a determinant of blood viscosity. Science 1970;168:977-9. |
| 9. | Easthope PL, Brooks DE. A comparison of rheological constitutive functions for whole human blood. Biorheology 1980;17:235-47. |
| 10. | Chien S. Biophysical behavior in suspensions. In: Surgenor D, editors. The Red Blood Cell. New York: Academic Press; 1975. |
| 11. | Bucherer C, Lacombe C, Leli'JC LL. Viscosity of the human blood. In: Jaffrin MY, Goubel F, editors. Fluid and Tissue Biomechanics. Paris: Masson; 1998. |
| 12. | Thiriet M, Martin-Borret G, Hecht F. Shear-thinning flow in a bend and a symmetrical planar bifurcation. Application to the blood flow in large vessels. J Phys III 1996;6:529-42. |
| 13. | Iolov A, Kane A, Bourgault Y, Owens R, Fortin A. A finite element method for a microstructure-based model of blood. Int J Numer Method Biomed Eng 2011;27:1321-49. |
| 14. | Thiriet M, Graham JM, Issa RI. A pulsatile developing flow in a bend. J Phys III 1992;2:995-1013. |
| 15. | Womersley JR. Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known. J Physiol 1955;127:553-63. |
| 16. | Quarteroni A, Ragni S, Veneziani A. Coupling between lumped and distributed models for blood flow problems. Comput Vis Sci 2001;4:111-24. |
| 17. | Boissonnat JD, Chaine R, Frey P, Malandain G, Salmon S, Saltel E, et al. From arteriographies to computational flow in saccular aneurisms: The INRIA experience. Med Image Anal 2005;9:133-43. |
| 18. | Vétel J, Garon A, Farinas MI, Thiriet M. Steady State Carotid Flow in vitro PIV measurements. 8 th World Congress on Computational Mechanics (WCCM8) and 5 th European Congress on Computational Methods in Applied Sciences and Engineering (Eccomas'08). Venice (Italy); 2008. Available from: congress2.cimne.com/eccomas/proceedings/eccomas2008. [Last accessed on 2019 Nov 12]. |
[Figure 1], [Figure 2], [Figure 3]
[Table 1], [Table 2], [Table 3]
|
Search Pubmed for