1. Introduction
Small diameter blood vessels followed by upstream vasoconstriction is often prone to secondary stenosis [
1]. Hence, regenerated blood vessels will be an appropriate solution for replacing these small diameter vessels (which were affected with secondary stenosis). Various techniques, such as extrusion, electrospinning, thermal-induced phase separation, braiding, hydrogel tubing, and gas foaming had been used for developing artificial blood vessels [
2,
3,
4,
5]. Serpentine vascular geometry having a rigid wall has been used in vascular regeneration [
6]. Unlike a straight channel blood vessel, both bifurcated and arch-type blood vessels do not facilitate constant and unidirectional blood flow [
7]. Instead, the flow becomes turbulent with larger eddies and torsion. Such phenomena encourages formulation and analysis of the flow driven mechanical force distributions within the serpentine vascular network.
Bioprinting methodologies commonly use hydrogels as a cell-laden material as they are able to come into contact with the cells without damaging their viability [
8,
9,
10,
11]. Moreover, hydrogels can readily be mixed with cells, while simultaneously allowing for high cell density and homogenous cell distribution throughout the scaffold [
12]. The source of the hydrogel is an important factor because of varying chemical and mechanical properties of synthetic and natural polymers [
13]. Naturally derived hydrogels, such as GelMA, consist of inherent signaling molecules that promote cell adhesion. On the other hand, hydrogels derived from other organisms such, as alginate, lack these binding sites for cell adhesion and attachment [
14]. Natural biomaterials regulate interactions between the cells and the extracellular matrix and provide an excellent environment for the cells to grow, differentiate, and proliferate [
13,
15]. Since GelMa is a gelatin based bioink, it can also resolve intricate vascular networks and channels that offer endothelial cells with the essential properties of their native environments. Rheological modifiers also called rheological additives are mixed with gelatin based bioinks to enhance the rheology of the resulting bioink, especially with regard to viscosity and yield stress.
The pulsatile and sinusoidal flows have been reported as the simplified form of physiological blood flow in circulatory system [
16,
17,
18]. These inlet velocity profiles represent simulated biomimetic flow conditions for the serpentine vessels [
19,
20,
21,
22]. Further, the flow facilitates the seeded endothelial cells with the physiological stress at the luminal wall of the serpentine structure [
23]. Limited study can be found dealing with both physiological and sinusoidal flow to generate hemodynamic stress on the annular surface of the straight and bifurcated substrates [
24,
25,
26]. However, such inlet boundary conditions were not considered in a serpentine network with an elastic wall. Previous studies had considered serpentine models independent of the upstream stenosis (or secondary stenosis). The physiological flow of blood had been found to reproduce adequate stress for the seeded endothelial cells on the luminal wall, while the rigid model of the wall compensates the viability of such cells under a shear stress condition [
27,
28]. This stress was generated due to the flow of culture media through the bioprinted channels. Due to the positioning of endothelial cells, it is more sensitive and responds quickly to the fluid flow driven stress [
29]. It had also been reported that generation of vascular growth factors, conservation of blood vessels and migration of endothelial cells were being improved in the presence of modulating wall shear stress [
30,
31]. Shear stress derived from flow structures were also found to influence interconnected points between the cytoskeleton and extracellular structures of endothelial cells [
32,
33,
34]. However, it was found that elasticity of the vessel wall promotes migratory phenomena of endothelial cells in the presence of lateral stress in the periodic form of lymphatics. Hence, it is very vital to consider the bioprinted serpentine structure with the elastic wall for enhancing the cell viability of seeded cells in the presence of physiological and sinusoidal flow conditions.
Most of the previous studies deal with the systemic development of a 3D scaffold for formulation of straight channel blood vessels. In addition, a perfusion bioreactor has traditionally been used because of its mechanical advantages rather than a rotatory bioreactor. The present study deals with the use of a 3D printer along with a rotatory bioreactor to fabricate a serpentine vascular structure. The use of a rotatory bioreactor facilitates adequate amounts of shear stress for stimulating the growth of human endothelial cells on the vessel walls. Numerical analysis was conducted to enumerate the relationship between flow driven different mechanical parameters and different inlet and wall boundary conditions. Variation of pressure, stress, and shear rate were found influenced by the inlet and outlet boundary conditions as well as by the viscosity of the acting fluid. There is limited data on cellular viability of vascular channels bioprinted using GelMA bioink mixed with PEGDA photoink. Cellular viability of GelMa hydrogels in a 3D cell culture model demonstrated that the hydrogel scaffolds provide a cell promoting environments for mesenchymal stem cells [
35]. Experimental characterization of the bioprinted channels has been performed to validate the optimum flow-driven mechanical parameters, while maintaining higher cellular viability.
3. Results
Numerical simulation was performed for evaluating flow induced mechanical parameters (stress, pressure, shear rate, and velocity) on the inner lumen of the serpentine model. For condition I (slip wall with sinusoidal flow, refer to
Table 3), significant variation of pressure, shear rate, velocity, and stress were observed across the entire model, as shown in
Figure 6. For velocity models V
1 and V
3, there was no variation in the pressure parameter over time. Maximum pressure generation for velocity models V
1, V
2 and V
3 were recorded as 2–5 × 10
5 Pa. Similarly, minimum pressure generation for different velocity models V
1, V
2, and V
3 were obtained as 0.2 × 10
5, −1 × 10
5, and −1 × 10
5 Pa, respectively. In the velocity model V
2, pressure generation on point P
1 at a time T was maximum, compared to other points, as shown in
Figure 6a. It has been observed that the shear rate profile was globally constant for all velocity models, ranging from 0.5 × 10
4 to 4.5 × 10
4 1/s. The shear rate results for all velocity models, around points P
1 and P
3 at time T was maximum, as shown in the
Figure 6b. For all velocity models, it can be concluded that the flow development was linearly correlated with the progression of time within the T cycle, as shown in
Figure 6c. The minimum and maximum axial velocity of 2 mm/min at neck region (P
1) and 20 mm/min in between the abdominal (P
2) and rear regions (P
3), respectively, were observed for condition I. From
Figure 6d, it had been inferred that the local maximum stress having a magnitude of 18 × 10
−10 N/m
2 had been generated for the case of velocity model V
2 at the end of the cycle in between the abdominal (P
2) and rear regions (P
3).
For the sinusoidal flow with no-slip wall boundary condition (condition II), local minimal variation in pressure, shear rate, and axial velocity were analyzed for velocity models V
1, V
2, and V
3 (
Figure 7). A significant variation of pressure distribution was recorded for velocity model V
1 compared to the two other velocity models. For velocity models V
2 and V
3, the local maximum pressure equal to 4 × 10
5 Pa was observed in the neck region (P
1) of the channel at end of cycle time T, as shown in
Figure 7a. Maximum pressure was generated at the time T/2 sec in the periphery of the rear region (P
2). It was noticed that the value of axial velocity was similar for V
1, V
2, and V
3 velocity models. The maximum velocity was generated at the end of the T cycle for all inlet velocity models. The minimum and maximum axial velocities were observed as 0 and 20 mm/min, respectively, in the neck region (P
1) and periphery of the rear region (P
3), as shown in
Figure 7c. Negligible variation had been observed in the shear rate for different velocity models in all regions of the serpentine at T/3, T/2, and T sec. The shear rate profiles for velocity models V
1, V
2, and V
3 were similar and its range of magnitude was 0.5–4.5 × 10
5 1/s, as shown in
Figure 7b.
The contour plot of pressure for velocity models V
4, V
5, and V
6 at a measured time period T/3, T/2, and T sec was evaluated similarly for physiological flow with no slip wall boundary condition (condition III). For each velocity model, it was observed that local maximum pressure (3.5–4 × 10
5 Pa) developed at the T/3 and T/2 sec at the neck region of the channel, as shown in
Figure 8a. For the same condition, the local maximum (at neck region) and minimum (rear and abdominal regions) shear rate value were predicted as 4.5 × 10
5 1/s and 0.5 × 10
5 1/s, respectively (
Figure 8b). The range of local minimum and maximum axial velocity was calculated between 0 to 18 mm/min, respectively (
Figure 8c). The velocity contour was found to be higher for periods T/3 and T/2 sec. Maximum velocity had been developed at the neck region of the serpentine channel.
The contour plot of physiological flow with slip wall boundary condition (condition IV) depicts that for velocity model V
4, the pressure magnitude at T/3 and T/2 sec was minimum (0.2 × 10
10 Pa). The maximum pressure (1 × 10
10 Pa) was developed at approximately the periphery of the neck region and the rear region at the end of the T cycle, as shown in
Figure 9a. At T/3 and T/2 secs, local maximum pressure was developed near the periphery of P
1, as shown in
Figure 9a. The contour plot of the shear rate was found similar for all velocity models of the physiological flow condition. The minimum shear rate corresponding to V
4, V
5, and V
6 were obtained as 1 × 10
7, 0.2 × 10
5, and 0 × 10
5, respectively. Similarly, the maximum shear rate for V
4, V
5, and V
6 were 8 × 10
7, 1.8 × 10
5, and 2.5 × 10
5, respectively, as shown in
Figure 9b. A maximum axial velocity of 200 mm/min was developed on all measured points at T sec, which was extremely high as compared to other conditions shown in
Figure 9c. For velocity models V
5 and V
6, a maximum pressure of 2.5 × 10
5 Pa had been observed on the entire channel at T/3 and T/2 sec at the inlet region. For the same period, the local maximum velocity for model V
6 were found to be higher than that of V
5 in the neck as well as the periphery of rear region. For the results obtained in
Figure 9d, it was observed that the stress decreases with an increase in velocity. Global values of stress for all velocities fluctuated between 2 and 25 × 10
−10 N/m
2 in the rear region of the serpentine channel.
After evaluating the magnitude of various mechanical parameters, the printed serpentine channels were imaged every alternate day using fluorescence microscopy to assess cell attachment and cell distribution within the channel. A sample channel image stained with CFSE is shown in
Figure 10a, but the CFSE fluorescent dye may still be producing fluorescence even if the cells are dead. So, for longer time periods samples were imaged without CFSE dye. Sample image of the channel and a cross-section of the channel after 14 days is shown in
Figure 10b. Analysis of the LIVE-DEAD assay demonstrated that on average 98.5% of the cells were viable after 14 days in the incubator, as shown in
Figure 10c, corresponding to the inlet flow boundary condition of 4.62 mm/min for the luminal fluid (or media).
Figure 10d, show a live dead assay fluorescence image for an inlet boundary condition of 4.48 mm/min and 4.76 mm/min, respectively. It is evident from
Figure 10d, that cell density in the serpentine channel is less than that corresponding to optimized mechanical parameters, as presented in
Figure 10c.
4. Discussion
The variation in magnitude of stress, pressure, shear rate, and axial velocity of the fluid flowing through the serpentine vascular channels depends on the material properties, orientation of serpentine geometry, and the flow architecture of working fluid at points P
1, P
2, and P
3. For condition I (as given in
Table 3), the proposed architecture follows the Hagen–Poiseuille model for a Newtonian fluid where pressure and velocities are proportional to each other [
44]. Variation in the axial velocity influences the shear rate distribution in a closed channel flow. The magnitude of the shear rate describes the flow behavior of working fluid inside the serpentine channel. A minimal variation in shear rate has been found due to a nominal variation in the axial velocity with respect to the radius of the serpentine because of considered rigid wall configuration [
45]. Variation in axial velocities also influences the downstream generated pressure in serpentine structures. Localized maximal pressure was generated due to chaotic motion of fluid particles. When the axial velocity of a fluid increases, some of the energy used by the random motion particles to follow the fluid direction develops a lower downstream pressure [
46]. The developing downstream pressure further gives rise to localized stress. It was also correlated that the velocity enhancement at the downstream region leads to formation of maximal stress. At higher axial velocity, fluid flow moved slowly near the wall due to diffusion and dispersion of fluid particles [
47]. Hence, the inlet velocity V
2 was found to create higher stress than its counterpart velocity V
3 = Performing transient analysis, the maximum variation in the axial velocity profile was obtained at the end of the full cycle in condition II (see
Table 3) due to positive acceleration of fluid particles [
48].
For the no-slip wall boundary condition with physiological flow at the inlet of serpentine (as given in
Table 3), a small variation compared to the maximum magnitude of the shear rate was obtained near the wall due to rigidity of the wall of the serpentine channel as well as due to negligible variation of axial velocity magnitude [
49]. There was minimum fluid velocity near the wall, which produces a non-significant variation of axial velocity near the wall. According to Hagen–Poiseuille model, maximum pressure magnitude was responsible for producing maximum axial velocity in a serpentine vascular model [
50]. For a no-slip condition, when fluid comes into contact with the wall of the serpentine channel, there is no relative movement between them, which decreases the stress magnitude [
51].
For the slip boundary wall condition with physiological flow at the inlet of the serpentine (as presented in
Table 3), the abrupt changes in axial velocity were contoured by the inlet velocity profile (V
4), which occurred due to slip conditions at the end of the T cycle. When a fluid flow comes into contact with the curvature section of the serpentine, it creates Coriolis force, which causes a transversal slope in the flowing fluid. Resultant interaction between Coriolis force and transversal slope develops a secondary force on the flow cross-section. This secondary force had disseminated to the curvature section, producing a higher axial velocity magnitude [
52]. Maximum axial velocity and presence of curvature of proposed serpentine structure were the effective parameters for developing maximum pressure near the neck and rear region of the serpentine structure at the end of T cycle.
The temporo-spatial deviations in flow-derived parameters were found to affect cell viability and functionality, as given in
Table 7. For condition I, the overall maximum and minimum deviation was obtained for pressure and shear rate parameters, respectively. The minimum deviation was obtained due to small changes in the measured point’s value. Serpentine structure and different axial velocity magnitude were responsible for producing maximum pressure deviation [
53]. The maximum velocity deviation was obtained due to development of Coriolis force by the serpentine structure, which had caused a radial pressure gradient [
54]. Increased pressure gradient radially had further influenced the radial flow of the media, thereby enhancing the deviation in shear stress of adjacent regions [
55].
For condition II (refer to
Table 3), minimum deviation was obtained for the axial velocity on point P
2 at the end of T/3 cycle. Similarly, the minimum deviation was reported for the pressure on point P
3 at the end of the T cycle. As the flow approached near the wall, it was converted into a transitional flow due to elasticity of the wall of the serpentine. Hence, the velocity of flowing fluid increases, producing a maximum deviation of velocities [
54].
The laminar flow used in the model and the absence of resistivity to the flow of fluid particles were responsible for producing a slight deviation in the velocity profile [
56,
57]. For condition III, it was observed that the deviation of all parameters had increased with respect to time. The variation in measured parameters had occurred due to the serpentine structure of the model, responsible for inducing the heterogeneity in the nature of the flow [
58]. The maximum and minimum deviation was obtained at point P
3 in boundary condition IV (refer to
Table 3) due to the physiological relevance flow in the presence of the elastic wall. Such a flow condition brought a non-stationary axial velocity profile over the period of time [
59].
The numerical analysis aided in the selection of optimized values of different biofluid dynamics features, such as axial velocity, shear rate, pressure, and stress of the proposed serpentine blood vessel model for producing an ideal condition for cell viability. High cellular viability was observed even after 14 days of seeding of HUVECs in the serpentine channel maintained in the bioreactor. Cell attachment was almost uniform across the channel as can be observed from the fluorescent microscopy and confocal images.