Shear Stress Quantification in Tissue Engineering Bioreactor Heart Valves: A Computational Approach

1. Introduction

Heart valve disease affects 2.5% of the US adult population, with higher incidence in older individuals. Currently, treatment options are limited and surgery is recommended for severe, symptomatic cases[1,2]. Valve replacement options include mechanical and bioprosthetic heart valves, but both have significant drawbacks. Mechanical heart valves are durable but require lifelong anti-coagulation therapies and result in non-physiological flow patterns. They have a higher risk of major adverse cardiovascular complications, particularly major bleeding[3]. Bioprosthetic heart valves can overcome some of the limitations of mechanical heart valves, but their durability is limited to 10-15 years, making them unsuitable for younger patients [4]. In pediatric patients, both types of valve prosthesis require frequent resizing surgeries due to their inability to remodel and adapt to growth.

Tissue engineering presents a promising opportunity to create heart valve prostheses that are superior to the existing options. A commonly used approach for developing tissue-engineered heart valves involves pre-seeding cells on a scaffold and maturing them in a bioreactor under controlled physiological conditions [5]. A properly matured tissue engineered heart valve (TEHV) seeded with human cells has the potential to grow, repair, and remodel once implanted. To achieve functional valve tissue, it is important to subject the seeded human cells to physiological conditions for maturation. There are different types of bioreactors that can produce in-vitro environments with different flowrate, pressure, temperature, and oxygen diffusion, to condition the TEHV[6]. Bioreactors can be broadly classified into flow-based, strain-based whole valve conditioning, and isolated cusp stimulation bioreactors[6,7]. The optimal maturation protocol for nurturing seeded cells on a valve scaffold is still unknown, and the mechanical cues resulting from flowrate and pressure differential are difficult to quantify. TEHVs are usually conditioned between 2-to-6-week cycles with mixed outcomes [6,8,9].

Previous attempts on valve conditioning have highlighted the importance of understanding the magnitude of wall shear stress (WSS) as the flows within the bioreactor during the conditioning protocol are progressively adjusted. VeDepo et al. investigated the potential for leaflet matrix restoration and repopulation using bone marrow-derived mononuclear cells; human aortic valves were decellularized and seeded with marrow-derived mononuclear cells and the test samples were subjected to different conditioning protocols. The study found that extended conditioning protocols (i.e., mechanical cues) led to retraction of valve leaflets, suggesting that an extended period of mechanical conditioning is not feasible [10]. The finding of VeDepo et al. suggest that the application of excessive WSS could potentially have a negative impact on the development and function of TEHVs. Kennamer et al. sought to investigate whether stem cells would differentiate into valvular interstitial cells when subjected to mechanical stimuli. The test group samples were subjected to mechanical cues for four weeks of in-vitro bioreactor conditioning, and the study found that most of the seeded cells had died, resulting in a large island of cell debris[11]. The study concluded by emphasizing the sensitivity of stem cells to mechanical cues and the need for careful, progressive adaptation to these cues.

The regulation of homeostasis in the endothelial cell lining depends significantly on WSS. Several ex-vivo studies [12,13,14] have shown that changes in WSS can lead to the valve becoming diseased. Although some studies [15,16,17,18,19] have provided an estimate of the WSS in the aortic valve during normal flow conditions, there is a lack of information abou the WSS distribution in TEHV and how it changes as flow increases during conditioning. Previous studies on heart valve WSS also assumed aortic valve geometries with sinuses (which may not be present in bioreactors) and considered blood-like properties to model fluid properties; however, culture media fluids in bioreactors have properties more similar to water. There is a need for obtaining WSS estimates in TEHV at various bioreactor flow rates. This knowledge would benefit researchers working on conditioning protocols, enabling them to optimize flow ramp up rates and gain a better understanding of cell adhesion and growth.

The objective of this study is to quantify WSS in TEHV during the maturation process in the bioreactor. To achieve this, we used fluid-structure interaction (FSI) simulations to generate valve geometries at different points in the pulsatile flow cycle; we then conducted computational fluid dynamics (CFD) simulations with these valve geometries to quantify shear stresses on the valve leaflets at multiple flow rates. The study focused on the aortic valve with annulus diameters of 12.3mm, 18.45mm, and 24.6mm, representing pediatric, adolescent, and adult heart valves, respectively.

2. Materials and Methods

The dynamics of heart valves are influenced by the fluid flow and pressure differences across them. FSI simulations can model valve dynamics by coupling the loading on the valve caused by the fluid to the disturbance on the fluid by valve motion. To achieve this, FSI simulations solve the continuum equations for the solid valve, as well as the continuity and Navier-Stokes equations for the fluid, to model the two-way fluid-structure interaction (between the valve and the fluid). Two approaches exist in modeling the FSI: 1) The monolithic approach which uses a single algorithm to solve the fluid and solid domain equations, 2) The partitioned approach which solves the two set of domain equations in separate algorithms and uses a coupling algorithm to balance the forces exchanged and continuity of the domain between the fluid and structural regions. This study uses the partitioned method [20,21] which has two distinct, overlapping meshes. We adopted the non-boundary fitted method where the fluid mesh remains fixed while the solid mesh deforms and intersects the fluid mesh, resulting in significantly reduced computational expense compared to traditional Arbitrary Lagrangian Eulerian methods where both meshes must deform. However, the limitation of this approach is that the fluid mesh being homogenous prohibits near-wall mesh refinement for obtaining accurate WSS quantification. We overcome this limitation of the FSI by extracting open valve geometries at different time points during systole, then using these geometries to perform subsequent CFD simulations which allow for boundary layer meshing to obtain detailed WSS. The details of each the FSI and the CFD phase are outlined in the sections to follow.

2.1. FSI Simulations

We use an idealized aortic valve with a uniform leaflet thickness of 0.4 mm, identical deformable leaflets, and a rigid corona with a ring at the base. The valve diameter and the ring external diameter is 24.6 mm and 28 mm, respectively; the total height of the geometry is 17 mm. The geometries for the smaller sized valve (12.30 and 18.45 mm) are obtained from scaling the reference model. Three idealized valves represent pediatric, adolescent, and adult sizes 12.3, 18.45 and 24.6 mm respectively (Figure 1). The non-boundary fitted FSI simulations we used is based on the study by Luraghi et al.[21] and has been validated experimentally for a polymeric heart valve in [22]. The study [22] demonstrated the accuracy of FSI simulation in estimating the valve kinematics such as the geometric orifice (GOA) against experimental results. The material model for the leaflet is linear elastic for the working strain range, with an elastic modulus of 3 MPa, Poisson’s ratio of 0.40, and a density of 1100 kg/m³. The fluid in the FSI simulation is modelled as a Newtonian fluid with a density of 1060 kg/m³ and a dynamic viscosity of 3 cP, similar to Luraghi et al.[21]. A physiological pressure difference is the driving boundary condition. We conducted the simulations over two cardiac cycles each with a duration of 0.8 seconds and extracted the results from the last period. All FSI simulations are performed on 8 cores of a 14 core Intel-MPI Xeon64 processor using the commercial finite element solver LS-Dyna Release 11 (ANSYS Inc, Canonsburg, PA, USA)

2.2. GOA and flow rates selection

We selected five GOAs, ranging from 45% to 100% of GOA (normalized to maximal values during the cycle), from each of the three FSI simulation for subsequent CFD simulation (Table 1). GOA is the anatomical area of the aortic valve orifice and is obtained from image planimetry [23]. Planform views (Figure 2) of the valve with dimensional references are imported into Fiji ImageJ [24] and measured for GOA. The GOA for the 15 geometries (5 GOA x 3 Valve sizes) selected for this study approximately range from 39 to 279 mm² (Figure 3).

Table 1. Summary of geometries selected from FSI simulation results to be used in CFD simulations.

Valve Diameter (mm)	GOA (mm2)	% Max GOA
24.60	279.4	100
	255.4	91.41
	218.4	78.17
	171.9	61.52
	129.5	46.35
18.45	160.9	100
	148.7	92.42
	118.9	73.90
	97.42	60.55
	76.79	47.73
12.30	71.58	100
	64.21	89.70
	54.34	75.92
	46.72	65.27
	39.73	55.50

Table 1. Summary of geometries selected from FSI simulation results to be used in CFD simulations.

Valve Diameter (mm)	GOA (mm2)	% Max GOA
24.60	279.4	100
	255.4	91.41
	218.4	78.17
	171.9	61.52
	129.5	46.35
18.45	160.9	100
	148.7	92.42
	118.9	73.90
	97.42	60.55
	76.79	47.73
12.30	71.58	100
	64.21	89.70
	54.34	75.92
	46.72	65.27
	39.73	55.50

Figure 2. Isometric and planform (top) views of GOA.

Figure 2. Isometric and planform (top) views of GOA.

Figure 3. Summary of the geometries included in the CFD simulation to quantify WSS. Numbers represent the GOA in mm².

Figure 3. Summary of the geometries included in the CFD simulation to quantify WSS. Numbers represent the GOA in mm².

For each geometry extracted from the FSI simulation, we performed multiple CFD simulations with a range of flow rates. The purpose of this is to collect data corresponding to different potential combinations of GOAs and flow rates, considering the progressive increases of flow rates in a typical TEHV conditioning protocol in the bioreactor. Prior studies [15,16,25,26] reported instantaneous peak flow rates ranging from 20 – 25 L/min for an average cardiac output of 5 L/min. We thus selected the following flow rates for the CFD simulations: 5, 11, 18, 24 and 30 L/min, which encompasses the wide range of flow rates that could occur in the bioreactor. Table 2 summarizes the flow and GOA combinations we include in this study.

2.3. CFD Simulations

The extracted valve geometries are processed in a pre-processor package Hypermesh (Altair Engineering, Inc., USA) to create the two-dimensional surface mesh for use in CFD. This surface mesh consists of the valve housed inside a cylinder of diameter 32.4 mm, analogous to how a TEHV is housed inside a bioreactor. Using the tip of the valve as a reference point, the inlet to the valve is situated 60 mm upstream, representing the typical short bioreactor inlet, and the outlet is located 1200 mm downstream to ensure numerical stability. The geometric model thus forms an enclosed fluidic volume with the valve situated within. The size of the domain is the same for all valves we simulated.

This surface mesh is imported into ANSYS FLUENT 19 R2 (ANSYS Inc., Canonsburg PA, USA) for discretizing the flow volume with three-dimensional elements and generating 12 boundary layers on the ventricular side of the valve. The volumetric mesh is composed of hexahedral elements at the center with transition into polyhedral elements near the surface (poly-hexcore). Grid independence study is carried out for the largest GOA of each valve (see supplemental materials).

In the CFD simulations, we model the fluid properties to resemble water at 37 °C with a density of 993.3 kg/m³ and dynamic viscosity of 6.914x10^-4 Pa.s. Prior studies [15,16],[25,26] in an attempt to quantify stresses under physiological conditions in-vivo, model the fluid to resemble blood (density 1040~1060 kg/m³, dynamic viscosity 3.5 cP). However, our viscometry measurement (using Cannon 9721-R59 viscometer) of the culture media indicates that the media closely resembles water (density 1008.5 kg/m³, dynamic viscosity 9.141x10^-4 Pa.s). Thus, we decided to model the fluid to resemble water at 37 °C. Constant magnitude velocity (plug type) is prescribed at the inlet. This is justified as most bioreactors have short inlet lengths. The plug type inlet velocities are derived from volumetric flow rate based on the continuity equation: Q = V * A, where Q is the flow rate in m³/s, V is the mean velocity in m/s and A is the cross-sectional area of the cylinder domain (824.48 mm²). At the outlet, an “outflow” boundary condition is prescribed, which is recommended when the downstream pressure and velocity are unknown [27]. The solution method uses the transient pressure-based solver with a laminar model with cell based least squares for the gradient spatial discretization. The pressure and momentum equations use second order and second order upwind discretization, respectively. The transient formulation is a bounded second order implicit scheme. The solution method a Pressure-Implicit with Splitting of Operators (PISO) form of pressure-velocity coupling. Relaxation factors for the solvers are set to default, while the target residuals of continuity and velocity are set to 1x10^-6. The time step size is based upon the equation Δt = L/(3V) where L & V are the characteristic length and velocity, respectively. The time step size is set to 0.01 seconds to standardize all the simulations. All simulations use the standard form of solution initialization based on the inlet velocity and have a simulated flow time of at least 20 seconds. The simulations are performed on a mixture of 24 and 28 core compute nodes on the Palmetto cluster, which is Clemson University’s primary high performance computing resource.

3. Results

3.1. 24.6 mm Valve

When the flow rate is low (between 5-11 L/min), the WSS histogram has a shape like an impulse function (Figure 4). This implies that the entire leaflet surface is subjected to similar levels of WSS. As the flow rate and GOA increase, the histogram tends to become more uniform in shape; this implies that different areas of the valve experience drastically different levels of WSS. Despite this, all the histograms are right-skewed, indicating that most of the leaflets experience low levels of WSS (less than 40 Pa). Previous research [28,29] has demonstrated that shear stresses greater than 40 Pa can harm the endothelial lining.

3.2. 18.45 mm Valve

As compared to Figure 4, the smaller valve size (GOA) and higher velocities lead to increased stresses in the 18.45 mm valve (Figure 5). The graph clearly shows that as GOA decreases and flow rate remains constant, WSS increases. Although most of the leaflets experience less than 40 Pa WSS, some instances at low GOA and high flow rate show WSS between 40 to 70 Pa, which can potentially damage cells [28,29].

3.3. 12.3 mm Valve

In Figure 6, it is evident that the 12.3 mm valve experiences considerably higher levels of WSS when compared to the other valves. At a flow rate of 5 L/min, WSS is mostly <40 Pa. However, at higher flow rates, we have consistently observed WSS exceeding 100 Pa, which has the potential to cause erosion of the endothelial cells [28,29].

3.4. WSS distribution within the leaflet

We typically saw three regions of WSS in the simulation results (Figure 7a). The tip of the leaflet experiences high magnitude WSS; and especially near the corners of the leaflets, the highest magnitudes are observed due to fluid acceleration resulting from converging geometry. At the center of the leaflet, we observe a low magnitude region of WSS. At the base of the leaflet, moderate to elevated WSS are observed due to fluid re-circulation. This pattern is consistent across almost all flow rates and geometries, except for high flow rates of 12.3 and 18.45 mm valve, where we observe increased stresses at the center due to flow separation (Figure 7b).

3.5. Regression models of shear stress

Here we present a regression model to estimate the 50^th (median) and 99^th percentile WSS based on flow rate and GOA. We examine the 50^th percentile instead of the mean because the right-skewed nature of the histogram can result in a large difference between the two, and the knowledge of the WSS threshold which categorizes each half of the cell population is potentially more useful for TEHV maturation protocol development. The 99^th percentile data excludes influence from any numerical artifacts along the perimeter, especially near the corners of the leaflets, where the small gap in leaflets does not allow for a sufficient number of boundary layers in the model.

The best fit regression lines (Figure 8 and Figure 9) of the transformed WSS (50^th and 99^th percentile) versus a function of flowrate and GOA clearly show that the regression models are independent of valve size (i.e. the results from all three valve sizes fit to the same regression model). Equations 2 and 3 below describe the regression models for estimating the 50^th (median) and the 99^th percentile WSS, respectively.

Ln(τ₅₀) = 1.39 ln(41258 * Q/GOA)−13.40

(1)

ln(τ₉₉) = 1.47 ln(41258 * Q/GOA) – 12.85

(2)

where τ₅₀ is 50^th percentile (median) and τ₉₉ is the 99^th percentile WSS (Pa) respectively, Q is flowrate in m³/s and GOA is the geometric orifice area in m².

We found that the data points located close together on figure 8 also exhibit very similar WSS histograms (Figure 10).

4. Discussion

This research provides useful data quantifying shear stresses in TEHVs in terms of the valve GOA and flow rate. CFD simulations were used to quantify the WSS pattern and magnitude for three sizes of TEHVs. We presented geometrical patterns and histograms of WSS, and a regression model to estimate the 50^th and 99^th percentile WSS based on the GOA and flow rate.

We identified three regions of WSS variation in the leaflet: high WSS at the tip of the leaflet and the corners of the leaflets, low to moderate WSS in the center, and moderate to elevated WSS at the base of the leaflets caused by recirculating fluid. These findings agree with previous report of regional variation in WSS along the length of the valve (Fig 5a in ref [19]). For the smaller valves at high flow rates, there can be moderate WSS at the center due to detachment of the boundary layer in that region. This knowledge may aid in the design and maturation of TEHVs. For example, one may explore a staged seeding protocol where cells are first seeded in the high WSS areas to begin the maturation process in the bioreactor at low flow; as the bioreactor flow rate ramps up, additional cells can be seeded in the low WSS area. This could provide a means to apply a more uniform WSS history to all of the cells in the TEHV.

Prior research [15] has reported WSS of 7.9 Pa on a 23 mm polyurethane valve at a location midway between the leaflet coaptation area and the center of the leaflet, under a flow rate of 22.5 L/min and fluid viscosity of 3.8 cP and density of 1.05 g/cm3. In our simulation of the 24.6 mm valve at the fully open position and subjected to a flow rate of 24 L/min, the WSS in the same region was approximately 8 Pa. Another prior investigation [16] reported WSS values of 7.1 Pa at the center for a valve with 24.8 mm diameter at flow rate of 18 L/min. The fluid used had a dynamic viscosity of 3.5 cP. In our study, the 24.6 mm valve subjected to an 18 L/min flow rate showed WSS of approximately 2.5 Pa at the center of the leaflet. The WSS values we found in our study are similar or slightly lower compared to previous reports. Our lower WSS values are likely due to the fact that we considered a fluid with properties resembling water and typical cell culturing media, where the previous studies used fluids with higher viscosity to represent blood. Theoretically, given a specific velocity gradient, WSS scales linearly with the dynamic viscosity. Our results suggest that, with a less viscous fluid, the velocity gradients in our simulation results are similar to those in previous studies well within one order of magnitude.

In this study, we investigated 66 unique combinations of GOA and flow rates across three different valve sizes. Our results reveal the role of GOA as the primary geometric parameter influencing WSS estimates. This conclusion is supported by both regression equations and the WSS distribution histogram presented in this study. Several key observations substantiate this assertion. Despite the distinct 3D geometries inherent in different-sized valves, the data from valves of varying sizes conform to the same regression line. Furthermore, the WSS regression equations have dependence on only two factors, GOA and flow rate, with GOA being the only geometric parameter. Finally, the WSS distribution histogram is invariant for scenarios corresponding to closely located points on Figure 8, suggesting a direct connection between histogram results and their positions on the Figure 8 regression plot. The sensitivity of our reported WSS results exclusively to GOA as a geometric parameter implies that detailed variations of the precise valve shape is not paramount. This insight is important in interpreting and applying our reported data. When using this data to estimate WSS in a bioengineered valve, the exact geometry of the valve need not be identical to those in our simulation for our findings to be relevant and applicable.

The regression models we developed can help tissue engineering researchers quantitatively estimate the WSS a valve is experiencing by measuring the instantaneous flow rate and GOA of the TEHV in a bioreactor, irrespective of the valve size. For example, consider a hypothetical conditioning protocol on a 24 mm valve, where the flow rate is ramped to physiological levels over a 14-day period. Suppose that we measure the peak instantaneous flow rate and the corresponding GOA during the maturation period (Table 3). The 50^th and 99^th percentile peak stresses can then be estimated via equations 1 and 2. In this hypothetical conditioning example, we see that on day 5, half of the cells in the TEHV experience peak WSS <3.39 Pa, and half of the cells experience peak WSS between 3.39 Pa and 13.61 Pa; this scenario corresponds to the location on Figure 8 with x and y value of 10.51 and 1.22, respectively — the closest data point to this location is data point 12 in Table S1, and thus the WSS histogram for the 18.45 mm valve at 62% GOA and 5 LPM (Figure 5) can be used as an estimate. The calculations from the regression equations together with the WSS histogram and WSS geometric patterns provided in this paper enable the tissue engineer to obtain a detailed estimate of WSS levels across the TEHV surface over the maturation time period.

There are limitations regarding the use of this research’s data. The study considered idealistic aortic valve models; while all aortic vales have similar overall 3-leaflet shape, there may be differences in the geometry between valves such as the ratio of the length and width of the leaflets. While the valve geometry has the dominate effects on WSS, the surrounding geometry of the bioreactor may have measurable effects; in this study we considered a standard cylindrical shape of fluid chamber to represent the bioreactor and did not explore varying bioreactor shapes. The data of this study is valid for fluids that resemble the properties of water and standard cell culturing media; if the density and viscosity of the culturing media differ significantly from those used in our simulations, only the qualitative results such as the WSS geometric patterns and overall shapes of WSS histograms may remain useful. Finally, a future experimental validation of WSS in TEHV under varying bioreactor flow conditions will be useful in corroborating with the results of this study and solidify the WSS estimates.

In summary, this study contributes to the field of TEHV by offering a set of data enabling tissue engineers to estimate WSS in bioreactor valves, a crucial factor for optimizing cell maturation. The dataset includes WSS distribution patterns and histograms, complemented by robust regression models; this information can help shift the current paradigm of largely trial-and-error approaches to more intentional, informed adjustments to bioreactor cell maturation protocols. Our regression models reveal that WSS predominantly depends on GOA rather than the overall valve diameter, highlighting the potential for the translation of maturation protocols across different valve sizes.