Chronic kidney failure is associated with significant burden on patient and payor alike. Renal transplant is superior to dialysis with respect to survival, cost, and overall burden of care; however, transplant is severely limited by scarcity of donor organs. Approximately 63% of prevalent patients with end-stage renal disease (ESRD) require hemodialysis as a life-extending therapy,1 even though transplantation offers a survival advantage over hemodialysis.2 This persistent unmet need for donor organs has stimulated interest in technology-based approaches for overcoming the scarcity problem in organ transplant.3–7
Our team is using a biohybrid approach, combining a high-efficiency, silicon-chip hemofilter as a glomerular analogue and a bioreactor of renal tubule cells to accomplish the metabolic, regulatory, and transport functions of the kidney.6 This paper describes our approach to an iterative design process that identifies and eliminates potentially thrombogenic flow fields in candidate devices.
First, computational fluid dynamics (CFD) was used to model the flow field through a hemofilter blood path and predict areas of thrombus formation. The CFD results were verified using magnetic resonance imaging (MRI) to measure pulsatile flow through the device in vitro to simulate physiologic conditions. After iterative in silico and in vitro improvements, prototype blood flow paths were machined from polycarbonate. Because the purpose of this study was to minimize flow-induced thrombogenicity, the prototypes were machined to the dimensions of the intact hemofilter but did not contain the silicon nanopore filters. The prototypes were implanted in large animals to test whether locations of in vivo thrombus formation agreed with CFD predictions. In this study, we used CFD to model unsteady flow in a blood conduit (device A), and noticed agreement between regions of low wall shear stress (WSS) and thrombus formation. The shape of the blood conduit was revised to minimize regions of low WSS (device B), and we observed clot-free flow in 30 day implant studies. To our knowledge, this is the first report of physiologically relevant blood flow simulations used to minimize thrombogenic potential of a hemofilter design for renal replacement therapy.
Blood Conduit Design and Manufacturing
Two primary design criteria were defined for the implanted hemofilter device. First, the cross-sectional shape of the blood path needed to transition from the round conduits of vasculature and vascular grafts to the rectangular duct flow of a parallel-plate hemofilter, and back to vasculature. Second, to facilitate surgical anastomosis to artery and vein of similar size, a hilum-like configuration was sought: inflow and outflow conduits antiparallel but within a few centimeters of each other. These conditions are then further constrained by the need to prevent thrombosis and avoid stagnant or recirculating flow. The initial design transformed the flow cross-section between 6 mm diameter pipes and a 7 mm wide rectangular duct in a simple “U” shape (device A, Figure 1). Wall shear stress is the stress imparted on the conduit wall by the tangentially flowing fluid. Given the relationship between WSS and thrombus formation,8–11 this parameter is used in design of blood-contacting devices. Consequently, the shape was revised in response to in silico and in vivo findings to include a helical inlet12 and a gradual outlet curve (device B, Figure 1) to reduce the size and duration of low WSS regions. All flow paths were designed in SolidWorks 2008 x64 Edition (Dassault Systèmes, Vélizy-Villacoublay, France), and physical prototypes were milled from biocompatible (ISO10993) LEXAN, Resin HPS6 (Sabic; Riyadh, Saudi Arabia). Pieces were vapor polished and thermally bonded to form the blood flow path (Hayes Manufacturing Services, Sunnyvale, CA).
Flow Measurement by MR
Device A was imaged using a 63 mm inner diameter quadrature radiofrequency (RF) volume coil inside a Varian 9.4T horizontal bore MRI scanner (Agilent Technologies, Santa Clara, CA). A computer-controlled Masterflex L/S pump driver with dual Easy-Load II pump heads (Cole-Parmer, Vernon Hills, IL) drove pulsatile blood flow through the device, using LabVIEW System Design Software (National Instruments, Austin, TX) to prescribe cyclic flow parameters over a period of 560 ms (Figure 2). A pneumatic pillow was used to trigger magnetic resonance (MR) image acquisition. To reduce the longitudinal magnetization relaxation time, a 2 ml dose of gadopentetate dimeglumine (Magnevist; Bayer HealthCare Pharmaceuticals, Inc., Wayne, NJ) was added to approximately 500 ml of citrated bovine blood (Lampire Biological Laboratories, Pipersville, PA) pumped through the device, improving the signal-to-noise ratio of the images.
Phase contrast MR (PC-MR) data were acquired in the through-plane direction using a velocity-encoded spoiled gradient echo sequence for “axial” slices (i.e., bulk flow direction at the inlet and outlet was normal to the slice). Imaging parameters included the following: repetition time = 700 ms; field of view = 38.4 mm × 9.6 mm; slice thickness = 1 mm; voxel size = 0.15 × 0.15 × 1 mm3. The velocity encoding (VENC) parameter was 150 cm/s. To correct for background phase effects due to field inhomogeneity, two sets of images were collected for each VENC direction, with the polarity of the bipolar VENC gradients reversed between acquisitions. After the trigger signal, a series of images collected for 16 cine time points described the pulsatile flow at a temporal resolution of 35 ms.
After data acquisition, images and 2D phase maps were reconstructed using MATLAB (The Mathworks; Natick, MA) and ITK-SNAP (www.itksnap.org),13 and velocity maps were calculated after phase subtraction of the two acquisitions. Velocity encoding aliasing was addressed using a phase unwrapping algorithm.14 Through-plane MR velocity data were imported into tec360 (Tecplot, Inc.; Bellevue, WA) for comparison with simulation data.
In Vivo Experiments
All animal experiments were approved by Vanderbilt University’s Institutional Animal Care and Use Committee. Devices were implanted in twelve 25–30 kg, mixed breed, female class A dogs: five device A and seven device B (multiple channel heights; Table 1). The 2–3 cm long externally supported polytetrafluoroethylene (PTFE) grafts were attached by end-to-side anastomoses from aorta to device inlet and from inferior vena cava to device outlet. After blood flow was established, the device was secured to the psoas muscle adjacent to the inferior pole of the left kidney. Intravenous heparin (100 units/kg) was administered during the operation, with an additional 1,000 units administered when creating the arterial anastomosis. Postoperatively, the dogs received either Lovenox (0.5 mg/kg; Sonofi, Bridgewater, NJ) once each day or antiplatelet therapy with 1.5 mg/kg acetylsalicylic acid (ASA) per day (Table 1). Blood flow through the grafts was assessed weekly by Doppler ultrasound. Devices were explanted if flow was not detected, or at the end of planned experiment (24–31 days; Table 1).
Computational Fluid Dynamics Simulations
The computational geometry and grid were constructed using ICEM-CFD (ANSYS, Inc.; Canonsberg, PA), and laminar simulations were conducted in Fluent v. 15 (ANSYS, Inc.). Entrance and exit lengths were appended to the device geometry to ensure fully developed flow. The computational volume was discretized using tetrahedral elements and three layers of thin, rectangular prism elements at the near wall region to improve boundary layer and WSS estimations.
Grid independence was confirmed when an approximately 50% increase in element number between grids produced a relative error of maximum WSS less than 5% and a maximum velocity relative error at the outlet less than 3%. For all simulations, blood was considered a Newtonian fluid, which is reasonable for the scale of these devices,15,16 with a density of 1,053 kg/m3 and an absolute viscosity of 0.00368 kg/ms.
Simulation of in vitro experiment.
For simulation of the in vitro experiment, inlet boundary conditions were based on MR-measured velocities obtained near the inlet of the device. Cine velocimetry measurements were interpolated in MATLAB to produce a volumetric waveform describing blood flow (Figure 2). The median Reynolds number was 260, and the Womersley number was 3.1. Three cycles were concatenated, allowing investigation of temporal convergence. The volumetric flow rate was imposed as the inlet boundary condition of the simulation; the outlet was prescribed as having zero diffusion flux; and the geometry was considered a rigid boundary. Visualization of results was performed using Fluent and Tec360 (Tecplot, Inc.).
Comparisons were made between CFD-calculated, through-plane velocity maps and MR-measured velocity data. The CFD data were extracted from the inlet and outlet faces of the computational geometry. The MR data were acquired from the slice closest to the inlet/outlet of the device and positioned inside the device. Comparisons were made at five time points in the cycle at which the MR and exported CFD data aligned temporally; these included data at the initial time point and peak flow through the device.
Simulation of in vivo experiment.
For device A, representative, ultrasound-measured centerline velocity data were obtained in a graft immediately proximal to an implanted device. Discrete points from the measured velocity waveform were interpolated in MATLAB. The median Reynolds number was 600, and the Womersley number was 3.1. This waveform was used to approximate the flow rate17 and was imposed as the transient inlet boundary condition. The outlet was prescribed as having zero diffusion flux, and the geometry was considered a rigid boundary. Because resistance of the device is expected to change with geometry revisions, for device B, an inlet pressure waveform measured in dogs18 defined the inlet boundary condition; the outlet pressure was defined as 5 mm Hg, and rigid boundaries were assumed. The median Reynolds number/Womersley number pairs were 770/0.8, 1,330/1.2, 1,680/1.6, and 1,870/2.6 for the 500, 750, 1,000, and 2,000 μm channel width devices, respectively. The pulsatile simulation was performed over three cardiac cycles to verify temporal convergence.
Simulated and Measured Flows In Vitro
We compared velocity maps to test agreement between simulations and experiment. Figures 3 and 4 show inlet and outlet velocity contours at approximately the same locations for the in vitro and in silico models. Both CFD and MR results indicate a parabolic inflow velocity profile over time and a skewed outlet velocity profile. The mean maximum axial velocity errors were 8.9% and 10.4% at the inlet and outlet planes, respectively. The velocity maps demonstrate qualitative agreement between CFD simulations and in vitro measurements.
Simulated Flows and Thrombus Formation In Vivo
Given the relationship between shear and thrombogenicity, simulated time-averaged WSS (TAWSS) maps were constructed. The CFD-simulated TAWSS demonstrates locally low WSS regions at the outer wall of the proximal curve and at the inner wall of the distal curve (Figure 5). Because of high WSS values during systolic acceleration, the TAWSS maps obscure regions of persistently low WSS. Mean WSS patterns averaged over the cardiac cycle excluding systolic acceleration demonstrate larger zones of low WSS (Figure 6). In device B, the zones of low TAWSS are smaller than those seen in device A, and the location of low TAWSS “foci” varies with channel height (Figure 7).
In device A, clots developed in three of five implants; however, in device B, no clots were revealed upon explant (Table 1). Clots in two of the device A prototypes were identified as residing or originating at sites corresponding to regions of low TAWSS and low mean WSS (Figures 5 and 6) identified in the CFD simulations. Figure 8 shows an extracted clot with a platelet-dominated portion formed in the inner wall region of the distal curve. The location corresponds to CFD-simulated regions of low WSS, particularly over the nonsystolic acceleration portions of the cycle (Figure 6). This platelet-dominated, pale region likely marks the region of thrombus initiation.
In device B, no thrombi were observed in the experiments (Figure 7). Notably, although the simulated low TAWSS loci changed with channel height in device B, no clots developed in vivo, indicating the revised design performs well for the range of channel heights.
Device hemodynamics influence coagulation, platelet activation, and eventual thrombosis (e.g., 8–11,19,20). The accurate prediction of flow fields is an essential step in design of blood-contacting devices.8–11,21–28 Zones of low WSS can accompany recirculation areas or separation regions, beyond which flow reattachment convects platelets toward the surface.11 The combination of slow flow and physical forcing of platelet-surface interaction is a powerful stimulus for clot formation.
Distributions of WSS are a function of geometry, blood flow rate, and pulsatility. Early investigations demonstrated the importance of studying hemodynamics in the context of physiologic conditions.29 Our implantable bioartificial kidney relies on cardiac perfusion for blood circulation. We hypothesized that simulating in vivo conditions when integrating CFD in the device design cycle would reduce thrombosis. The simulations based on physiologic flows agreed with in vivo observed thrombus location. Design modifications eliminating areas of low WSS also eliminated device thrombosis over a 30 day implant period. To our knowledge, this is the first investigation of hemodynamics in a cardiac pressure–driven hemofilter using physiologic flow conditions.
Shear loading experienced by platelets traversing intravascular devices causes long-term platelet activation.30 Although the fluid shear stresses in our device are below hemolysis thresholds,6 subhemolytic levels of shear may prime platelets for local or distant thrombosis. Lagrangian and Eulerian approaches for calculating platelet damage28,31,32 and platelet residence time calculations may improve future thrombogenicity assessments.33
Our study is limited by the techniques used. Because of constraints of MR measurement and registration with CFD-calculated values, the locations of the inlet and outlet comparison planes were not identical in the measured and simulated data sets. Additionally, use of human viscosity values to simulate flow of citrated bovine blood may have resulted in velocity estimation errors. In the absence of in vivo pressure data, use of the diffusion flux boundary condition also may have contributed to discrepancies between measured and simulated velocity fields. A laminar model was used for these simulations. The vast majority of physiologic flows are laminar, but notable exceptions include arteriovenous grafts,34 severe stenosis,35 valve disease,36 and mechanical heart valves.37 Although instantaneous peak flows may reach Reynolds numbers associated with transitional flow under steady state conditions, they do not necessarily produce instabilities in pulsatile flow. Trip et al.38 (2012) demonstrated that for oscillating flow in a pipe, laminar flow persisted at mean Reynolds numbers higher than those investigated in this study. Further, they found that in the transition region, transitional flow was independent of Womersley number, contrary to previous works (e.g., 39,40). It should be noted that the Trip et al.38 study examined pipe flow. Although the geometries we considered are more complicated than pipe flow, our devices do not have sudden expansions.
The simulations presented here estimate hemodynamic parameters at high spatiotemporal resolution. Inclusion of in vivo and in vitro experiments allows demonstration of phenomena that models cannot predict, such as the accumulation of blood components on membranes and how this may induce a sieving effect or formation of clotting precursors. This underscores the utility of a combined in silico and in vivo/in vitro experimental approach when designing intravascular devices.
This study employed a combination of in vivo, in vitro, and in silico approaches to identify potential fluid dynamics contributors to thrombogenicity in an implantable hemofilter design under realistic operating conditions and to guide future device iterations. The results demonstrate both the feasibility and the utility of such a combined approach in developing intravascular, bioartificial, organ-replacement devices.
1. U.S. Renal Data System, USRDS 2016 Annual Data Report: Epidemiology of Kidney Disease in the United States. 2016.Bethesda, MD,
2. Wolfe RA, Ashby VB, Milford EL, et al.Comparison of mortality in all patients on dialysis, patients on dialysis awaiting transplantation, and recipients of a first cadaveric transplant. N Engl J Med 1999.341: 17251730,
3. Song JJ, Guyette JP, Gilpin SE, Gonzalez G, Vacanti JP, Ott HCRegeneration and experimental orthotopic transplantation of a bioengineered kidney. Nat Med 2013.19: 646651,
4. Abolbashari M, Agcaoili SM, Lee MK, et al.Repopulation of porcine kidney scaffold using porcine primary renal cells. Acta Biomater 2016.29: 5261,
5. Lanza RP, Chung HY, Yoo JJ, et al.Generation of histocompatible tissues using nuclear transplantation. Nat Biotechnol 2002.20: 689696,
6. Kensinger C, Karp S, Kant R, et al.First implantation of silicon nanopore membrane hemofilters. ASAIO J 2016.62: 491495,
7. Fissell WH, Manley S, Westover A, Humes HD, Fleischman AJ, Roy SDifferentiated growth of human renal tubule cells on thin-film and nanostructured materials. ASAIO J 2006.52: 221227,
8. Hochareon P, Manning KB, Fontaine AA, Tarbell JM, Deutsch SCorrelation of in vivo
clot deposition with the flow characteristics in the 50 cc penn state artificial heart: A preliminary study. ASAIO J 2004.50: 537542,
9. Topper SR, Navitsky MA, Medvitz RB, et al.The use of fluid mechanics to predict regions of microscopic thrombus formation in pulsatile VADs. Cardiovasc Eng Technol 2014.5: 5469,
10. Yin W, Alemu Y, Affeld K, Jesty J, Bluestein DFlow-induced platelet activation in bileaflet and monoleaflet mechanical heart valves. Ann Biomed Eng 2004.32: 10581066,
11. Duraiswamy N, Cesar JM, Schoephoerster RT, Moore JE JrEffects of stent geometry on local flow dynamics and resulting platelet deposition in an in vitro
model. Biorheology 2008.45: 547561,
12. Caro CG, Seneviratne A, Heraty KB, et al.Intimal hyperplasia following implantation of helical-centreline and straight-centreline stents in common carotid arteries in healthy pigs: Influence of intraluminal flow. J R Soc Interface 2013.10: 20130578,
13. Yushkevich PA, Piven J, Hazlett HC, et al.User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability. Neuroimage 2006.31: 11161128,
14. Ghiglia DC, Pritt MDTwo-Dimensional Phase Unwrapping: Theory, Algorithms, and Software. 1998.Wiley,
15. Pries AR, Neuhaus D, Gaehtgens PBlood viscosity in tube flow: Dependence on diameter and hematocrit. Am J Physiol 1992.263(6 pt 2): H1770H1778,
16. Lee SW, Steinman DAOn the relative importance of rheology for image-based CFD models of the carotid bifurcation. J Biomech Eng 2007.129: 273278,
17. Womersley JRMethod for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known. J Physiol 1955.127: 553563,
18. Olson RMAortic blood pressure and velocity as a function of time and position. J Appl Physiol 1968.24: 563569,
19. Moake JL, Rudy CK, Troll JH, et al.von Willebrand factor abnormalities and endothelial cell perturbation in a patient with acute thrombotic thrombocytopenic purpura. Am J Med Sci 1986.291: 4750,
20. Sakariassen KS, Joss R, Muggli R, et al.Collagen type III induced ex vivo thrombogenesis in humans. Role of platelets and leukocytes in deposition of fibrin. Arteriosclerosis 1990.10: 276284,
21. Lee JC, Lee K, Kim HCMathematical analysis for internal filtration of convection-enhanced high-flux hemodialyzer. Comput Methods Programs Biomed 2012.108: 6879,
22. Pelosi A, Sheriff J, Stevanella M, Fiore GB, Bluestein D, Redaelli AComputational evaluation of the thrombogenic potential of a hollow-fiber oxygenator with integrated heat exchanger during extracorporeal circulation. Biomech Model Mechanobiol 2014.13: 349361,
23. Rambod E, Beizai M, Rosenfeld MAn experimental and numerical study of the flow and mass transfer in a model of the wearable artificial kidney dialyzer. Biomed Eng Online 2010.9: 21,
24. Ge L, Leo HL, Sotiropoulos F, Yoganathan APFlow in a mechanical bileaflet heart valve at laminar and near-peak systole flow rates: CFD simulations and experiments. J Biomech Eng 2005.127: 782797,
25. Chiu WC, Slepian MJ, Bluestein DThrombus formation patterns in the HeartMate II ventricular assist device: Clinical observations can be predicted by numerical simulations. ASAIO J 2014.60: 237240,
26. Li Z, Kleinstreuer CAnalysis of biomechanical factors affecting stent-graft migration in an abdominal aortic aneurysm model J Biomech 2006.39: 22642273,
27. Loth F, Jones SA, Giddens DP, Bassiouny HS, Glagov S, Zarins CKMeasurements of velocity and wall shear stress inside a PTFE vascular graft model under steady flow conditions. J Biomech Eng 1997.119: 187194,
28. Morbiducci U, Ponzini R, Nobili M, et al.Blood damage safety of prosthetic heart valves. Shear-induced platelet activation and local flow dynamics: A fluid–structure interaction approach J Biomech 200942: 19521960,
29. Ku DN, Giddens DPPulsatile flow in a model carotid bifurcation. Arteriosclerosis 1983.3: 3139,
30. Bluestein D, Niu L, Schoephoerster RT, Dewanjee MKFluid mechanics of arterial stenosis: Relationship to the development of mural thrombus. Ann Biomed Eng 1997.25: 344356,
31. Marom G, Bluestein DLagrangian methods for blood damage estimation in cardiovascular devices—How numerical implementation affects the results. Expert Rev Med Devices 2016.13: 113122,
32. Hansen KB, Arzani A, Shadden SCMechanical platelet activation potential in abdominal aortic aneurysms. J Biomech Eng 2015.137: 041005,
33. Esmaily-Moghadam M, Hsia TY, Marsden ALA non-discrete method for computation of residence time in fluid mechanics simulations. Phys Fluids (1994) 2013.25: 110802,
34. Loth F, Fischer PF, Bassiouny HSBlood flow in end-to-side anastomoses. Annu Rev Fluid Mech 2008.40: 367393,
35. Lee SE, Lee SW, Fischer PF, Bassiouny HS, Loth FDirect numerical simulation of transitional flow in a stenosed carotid bifurcation. J Biomech 2008.41: 25512561,
36. Stein PD, Sabbah HNTurbulent blood flow in the ascending aorta of humans with normal and diseased aortic valves. Circ Res 1976.39: 5865,
37. Yoganathan AP, Chandran KB, Sotiropoulos FFlow in prosthetic heart valves: State-of-the-art and future directions. Ann Biomed Eng 2005.33: 16891694,
38. Trip R, Kuik DJ, Westerweel J, Poelma CAn experimental study of transitional pulsatile pipe flow. Phys Fluids 2012.24: 14103,
39. Nerem RM, Seed WAAn in vivo
study of aortic flow disturbances. Cardiovasc Res 1972.6: 114,
40. Falsetti HL, Carroll RJ, Swope RD, Chen CJTurbulent blood flow in the ascending aorta of dogs. Cardiovasc Res 1983.17: 427436,