Numerical and Experimental Analysis of Shear Stress Influence on Cellular Viability in Serpentine Vascular Channels

3D bioprinting has emerged as a tool for developing in vitro tissue models for studying disease progression and drug development. The objective of the current study was to evaluate the influence of flow driven shear stress on the viability of cultured cells inside the luminal wall of a serpentine network. Fluid–structure interaction was modeled using COMSOL Multiphysics for representing the elasticity of the serpentine wall. Experimental analysis of the serpentine model was performed on the basis of a desirable inlet flow boundary condition for which the most homogeneously distributed wall shear stress had been obtained from numerical study. A blend of Gelatin-methacryloyl (GelMA) and PEGDA200 PhotoInk was used as a bioink for printing the serpentine network, while facilitating cell growth within the pores of the gelatin substrate. Human umbilical vein endothelial cells were seeded into the channels of the network to simulate the blood vessels. A Live-Dead assay was performed over a period of 14 days to observe the cellular viability in the printed vascular channels. It was observed that cell viability increases when the seeded cells were exposed to the evenly distributed shear stresses at an input flow rate of 4.62 mm/min of the culture media, similar to that predicted in the numerical model with the same inlet boundary condition. It leads to recruitment of a large number of focal adhesion point nodes on cellular membrane, emphasizing the influence of such phenomena on promoting cellular morphologies.


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.

Methodology
A 3D model of the serpentine blood vessel was implemented using COMSOL Multiphysics. Results of the numerical models were further used to optimize the media flow parameters: pressure, shear rate, stress, and axial velocity for different inlet flow boundary conditions. These optimized parameters were further utilized for maintaining the printed serpentine channels in a rotary bioreactor. Finally, a microfluidic based chip system was developed to fabricate the serpentine vascular channels.

Numerical Analysis
The serpentine channel model is shown in Figure 1. Fluid flow (with Newtonian viscosity model for blood equivalent) was evaluated for different wall boundary conditions [slip (S) and no-slip (NS)]. To ensure that the vascular serpentine channel model reflects the realistic simulation of in vivo conditions, appropriate relevant boundary conditions were used for numerical analysis. The wall of the serpentine vascular channel was taken as a rigid body, with a no-slip boundary condition [36]. Fluid-structure interaction (FSI) physics was used for simulating the elastic behavior of the serpentine channel model in another study. For the simulation, the wall of the serpentine vascular structure was considered as a moving mesh. The following properties of blood were used for numerical analysis: density = 1066 kg/m 3 , dynamic viscosity = 3.5 × 10 −3 Pa.s. Young's Modulus = 0.4 × 10 6 Pa, Poisson's ratio = 0.5. GelMA material properties (density = 1020 kg/m 3 , dynamic viscosity = 4.2 × 10 −3 Pa.s, Young Modulus 3.18 KPa) were used for the modeling of the serpentine wall. Simulations were performed at a reference temperature of 293.15 K. Navier-Stokes equation was used to represent the momentum of the fluid model. Time-dependent partial differential equation of incompressible Navier-Stokes (Equation (1)) was used for the simulation.
where u is the blood velocity and P is the blood pressure, ρ is the blood density and is the blood viscosity. Fluid was considered to be incompressible and hence a continuity equation of incompressible form was used. Linear elastic material properties of solid mechanics and fluid properties of laminar flow were considered for the study. An isotropic solid model was used for both slip and no-slip conditions. A fully coupled option was selected in COMSOL Multiphysics to establish a possible connection between solid mechanics physics and laminar flow physics. In FSI analysis, all of the meshes created were physics controlled and automatically generated at a normal element size. The statistics of mesh for slip and no slip boundary conditions are listed in Table 1.
A transient model was used to realize the transport phenomena of the fluid in a 2D serpentine channel. A MUMPS based direct linear solver was used to solve the model with a damping factor of unity and a recovery damping factor on 0.75. Velocities corresponding to mean arterial pressure: 80 mmHg, 85 mmHg, and 90 mmHg were considered for the study in Microchannel [37][38][39][40].   The physiological functions (waveform) were derived from 4D-Laser Doppler data of blood flow across a cross-section of sub-clavicular artery of healthy subject. The study was approved by Institute Ethical Committee (NITRR/IEC/2021/12). Obtained temperospatial functions were interpolated thereafter using Curve-Fitting MATLAB toolbox version 2018a, Mathworks, India. Thereafter, the interpolated function with lowest RMSE value was selected for the study. The functions corresponding to inlet velocity models for V 1 to V 6 are given in Table 2. One complete cycle of both sinusoidal and physiological flow corresponds to one complete cardiac cycle (represented by T). Different study conditions were formulated based on the above boundary conditions: no-slip (wall) and sinusoidal inlet flow, slip (wall) and sinusoidal inlet flow, no-slip (wall) and physiological inlet flow, slip (wall) and physiological inlet flow (as tabulated in Table 3). A no-slip boundary condition will exclusively provide the influence of inlet flow type in the development of flow physics of media inside the serpentine. On the other hand, slip wall boundary condition will simulate the elasticity of the serpentine wall (similar to the natural blood vessel wall). Table 2. Velocity function for V 1 to V 6 .

Inlet Velocity Model Parameters Value Function
Sinusoidal Flow (SF) Physiological Flow (PF) Micromachines 2022, 13, x FOR PEER REVIEW 4 of 20 The physiological functions (waveform) were derived from 4D-Laser Doppler data of blood flow across a cross-section of sub-clavicular artery of healthy subject. The study was approved by Institute Ethical Committee (NITRR/IEC/2021/12). Obtained tempero-spatial functions were interpolated thereafter using Curve-Fitting MATLAB toolbox version 2018a, Mathworks, India. Thereafter, the interpolated function with lowest RMSE value was selected for the study. The functions corresponding to inlet velocity models for V1 to V6 are given in Table 2. One complete cycle of both sinusoidal and physiological flow corresponds to one complete cardiac cycle (represented by T). Different study conditions were formulated based on the above boundary conditions: no-slip (wall) and sinusoidal inlet flow, slip (wall) and sinusoidal inlet flow, no-slip (wall) and physiological inlet flow, slip (wall) and physiological inlet flow (as tabulated in Table 3). A no-slip boundary condition will exclusively provide the influence of inlet flow type in the development of flow physics of media inside the serpentine. On the other hand, slip wall boundary condition will simulate the elasticity of the serpentine wall (similar to the natural blood vessel wall).
Three different locations corresponding to positions P 1 (4, −2.5), P 2 (5.8, 2.4) and P 3 (7.7, −2.5) on the neck, abdominal, and rear region, respectively, were identified on the serpentine channel. These three regions are the downstream region immediately after curvilinear orientation of ascending and descending phases of the serpentine structure, which are more prone to Coriolis forces within the serpentine model. Flow parameters (pressure, shear rate, stress) over the entire cardiac cycle were evaluated at an interval of T/3, T/2, and T cycles.

Grid Convergence Test
A grid convergence test was conducted for optimizing the grid size of the flow domain. Grid Convergence Index (GCI) provides a uniform measure of convergence for grid refinement study [41]. The discretization of sinusoidal flow for velocity V 1 with no slip condition had been selected for three different element sizes. For all element sizes-pressure, shear rate and velocity at probe points P 1 and P 2 (as marked in Figure 1) were measured for the entire time cycle. In the current study, the grid refinement ratio r, which is equivalent to mean refinement ratio r mean of r 12 and r 23 were calculated using Equation (3a) where, r 12 = element value of Fine discretization element value of Normal discretization r 23 = element value of Normal discretization element value of coarse discretization Richardson extrapolation [42] introducing the p-th order method, as shown by Equation (3b): where f 1 , f 2 and f 3 are magnitude of derived parameters (e.g., velocity) at point P 1 and P 2 at different grid resolution at time T.
With the different parameters listed in Table 4, values of r and p were evaluated from Equation (3a,b), respectively, and the calculated grid convergence magnitude for two probe points is listed in Table 5. Further, the Grid Convergence Index (GCI) was evaluated using the Richardson extrapolation method based on estimated fractional error [41], as shown in Equation (4a): The safety factor (Fs) selected for this study was considered to be 1.25. It was observed that GCI for the fine grid (GCI 12 ) is relatively low if compared to the coarse grid (GCI 23 ), which demonstrates that the numerical simulation dependency on the cell size was minimum with a decrease in grid size. Figure 2 shows the plot between extrapolated value and the Richardson extrapolation estimation for centerline velocity. The y-axis of Figure 2 represents Micromachines 2022, 13, x FOR PEER REVIEW 6 of 20 equivalent to mean refinement ratio % &'() of % and % were calculated using Equation (3a) where, % = * element value of Fine discretization element value of Normal discretization % = * element value of Normal discretization element value of coarse discretization Richardson extrapolation [42] introducing the <-th order method, as shown by Equation (3b): Where , and are magnitude of derived parameters (e.g., velocity) at point P1 and P2 at different grid resolution at time T.
With the different parameters listed in Table 4, values of r and p were evaluated from Equation (3a,b), respectively, and the calculated grid convergence magnitude for two probe points is listed in Table 5. Further, the Grid Convergence Index (@AB) was evaluated using the Richardson extrapolation method based on estimated fractional error [41], as shown in Equation (4a): The safety factor (Fs) selected for this study was considered to be 1.25. It was observed that GCI for the fine grid (@AB ) is relatively low if compared to the coarse grid (@AB ), which demonstrates that the numerical simulation dependency on the cell size was minimum with a decrease in grid size. Figure 2 shows the plot between extrapolated value and the Richardson extrapolation estimation for centerline velocity. The y-axis of Figure 2 represents ƹ KL ,K , which is the difference in the values of a particular parameter at two different grid resolution and given by: It was observed that the extrapolated values for the velocity at point P1 and P2 did not change significantly on further decreasing the mesh size to represent finer grid resolution.
where, % = * element value of Fine discretization element value of Normal discretization % = * element value of Normal discretization element value of coarse discretization Richardson extrapolation [42] introducing the <-th order method, as shown by Equation (3b): Where , and are magnitude of derived parameters (e.g., velocity) at point P1 and P2 at different grid resolution at time T.
With the different parameters listed in Table 4, values of r and p were evaluated from Equation (3a,b), respectively, and the calculated grid convergence magnitude for two probe points is listed in Table 5. Further, the Grid Convergence Index (@AB) was evaluated using the Richardson extrapolation method based on estimated fractional error [41], as shown in Equation (4a): The safety factor (Fs) selected for this study was considered to be 1.25. It was observed that GCI for the fine grid (@AB ) is relatively low if compared to the coarse grid (@AB ), which demonstrates that the numerical simulation dependency on the cell size was minimum with a decrease in grid size. Figure 2 shows the plot between extrapolated value and the Richardson extrapolation estimation for centerline velocity. The y-axis of Figure 2 represents ƹ KL ,K , which is the difference in the values of a particular parameter at two different grid resolution and given by: It was observed that the extrapolated values for the velocity at point P1 and P2 did not change significantly on further decreasing the mesh size to represent finer grid resolution.

Experimental Analysis
A co-axial nozzle configuration of a 3D printer was used to print the serpentine channel. The printing bed had been located inside closed chamber housing, as shown in Figure  3. The construct was printed using a blend of GelMA photoInk purchased from CELLINK, Boston USA and PEGDA200 photoInk (purchased from CELLINK, Boston, MA, USA). This photo-ink blend equipped the final construct with cell attachment sites, high shape fidelity, and minimal construct swelling. CaCl2 was used as a chemical cross-linker during printing process of the serpentine. The printing environment parameters are presented in Table 6. It was observed that the extrapolated values for the velocity at point P 1 and P 2 did not change significantly on further decreasing the mesh size to represent finer grid resolution.

Experimental Analysis
A co-axial nozzle configuration of a 3D printer was used to print the serpentine channel. The printing bed had been located inside closed chamber housing, as shown in Figure 3.
The construct was printed using a blend of GelMA photoInk purchased from CELLINK, Boston USA and PEGDA200 photoInk (purchased from CELLINK, Boston, MA, USA). This photo-ink blend equipped the final construct with cell attachment sites, high shape fidelity, and minimal construct swelling. CaCl 2 was used as a chemical cross-linker during printing process of the serpentine. The printing environment parameters are presented in Table 6.

Endothelization of Vascular Network
A total of 50-500 µL of HUVEC suspensions (10 7 cells per mL) was injected via a micro-pipette to fill the channel. Simultaneously, the inlet and the outlet of the construct was sealed with a pinch-clamp. Channels were injected with HUVECs stained with CellTrace CFSE dye (ThermoFisher, Waltham, MA, USA). The device had been centrifuged at 40 rcf for 1 min with a slow acceleration and deceleration to let the cells settle to the bottom. Incubation at 37 °C facilitated the adhesion of cells to form the innermost layer of the vessels. After incubation for 30 min, the construct was placed over 180° for cell adhesion on the other side of the vessel. Finally, the cells were incubated for 5 h at 37 °C on a shaker.
A blood-equivalent working fluid was put in a 10 mL BD Plastipak syringe. The filled syringe was loaded in the programable Shenchen syringe pump and was operated in infusion mode. A syringe pump was connected to microfluidic accessories (one-way luer lock check valve, microfluidic fitting female luer lock adopter kit, barbed to male adopter and syringe tube) in a series to complete the flow driven system. The other end of the tube was connected to a fabricated microfluidic channel placed in the petri dish via a barbed male adopter, microfluidic fitting female luer lock kit, male luer lock kit, fitting adopter

Endothelization of Vascular Network
A total of 50-500 µL of HUVEC suspensions (10 7 cells per mL) was injected via a micropipette to fill the channel. Simultaneously, the inlet and the outlet of the construct was sealed with a pinch-clamp. Channels were injected with HUVECs stained with CellTrace CFSE dye (ThermoFisher, Waltham, MA, USA). The device had been centrifuged at 40 rcf for 1 min with a slow acceleration and deceleration to let the cells settle to the bottom. Incubation at 37 • C facilitated the adhesion of cells to form the innermost layer of the vessels. After incubation for 30 min, the construct was placed over 180 • for cell adhesion on the other side of the vessel. Finally, the cells were incubated for 5 h at 37 • C on a shaker.
A blood-equivalent working fluid was put in a 10 mL BD Plastipak syringe. The filled syringe was loaded in the programable Shenchen syringe pump and was operated in infusion mode. A syringe pump was connected to microfluidic accessories (one-way luer lock check valve, microfluidic fitting female luer lock adopter kit, barbed to male adopter and syringe tube) in a series to complete the flow driven system. The other end of the tube was connected to a fabricated microfluidic channel placed in the petri dish via a barbed male adopter, microfluidic fitting female luer lock kit, male luer lock kit, fitting adopter male luer to male thread, and mini luer to luer adopter. All the microfluidic components were procured from DARWIN microfluidics.
HUVECS media was fed to the reservoir for cell culture over the printed layers. Regulatory control of each component, such as media pressure, temperature, valve, pH and CO 2 and N 2 delivery, etc., were achieved. Induction of media to the housing was made by using a pulsatile pump, as shown in Figure 4. It creates pulsatile hydrostatic pressure to the printed bioink construct, similar to the optimized inlet velocity profiles finalized from numerical study. The printed constructs were imaged using fluorescence microscopy to evaluate cell attachment and cell distribution within the channel.
Micromachines 2022, 13, x FOR PEER REVIEW 9 of 20 male luer to male thread, and mini luer to luer adopter. All the microfluidic components were procured from DARWIN microfluidics. HUVECS media was fed to the reservoir for cell culture over the printed layers. Regulatory control of each component, such as media pressure, temperature, valve, pH and CO2 and N2 delivery, etc., were achieved. Induction of media to the housing was made by using a pulsatile pump, as shown in Figure 4. It creates pulsatile hydrostatic pressure to the printed bioink construct, similar to the optimized inlet velocity profiles finalized from numerical study. The printed constructs were imaged using fluorescence microscopy to evaluate cell attachment and cell distribution within the channel.

In-Vitro Testing
Microfluidic channel had been formulated by combining two casted segments. The upper part and the lower part of the printed structure of the model were merged to form a microfluidic system to support a serpentine structure, for housing the printed construct. Media had been transfused through the serpentine structure from the inlet of the structure in the form of sinusoidal waveform and physiological waveform using a programmable syringe pump (DARWIN microfluidic Shenchen Syringe pump, Paris, France), as shown in Figure 4.

Uncertainty Analysis
Uncertainty analysis was performed to validate the model characteristics and to determine the shear stress sensitivity. The uncertainty equation of shear stress on the fabricated channel model was represented by Equation (5). Shear stress was generated when the printed serpentine model came into contact with the flowing HUVECS media.
Here, W X , W ] and W [ are the standard deviations of the pressure, velocity, and the strain rate standard deviation, respectively. From Equation (5), the uncertainty value of stress equivalent was evaluated as ±179.36501 N/m 2 .

Cellular Viability Assessment
The viability of the bioprinted vascular channel was evaluated over a period of 14 days by performing LIVE-DEAD assay. The channels were stained with the LIVE-DEAD cell imaging kit (Fisher Scientific, Waltham, MA, USA). Live cells were stained with Calcein-AM (2 µM), the fluorescence was measured at 494-517 nm, and dead cells were stained with ethidium homodimer-1 (4 µM), which was measured at 528-617 nm. The channels were incubated with LIVE-DEAD stain for 30 min and the samples were washed with PBS. Imaging was then performed using a Nikon C1 Confocal fluorescence

In-Vitro Testing
Microfluidic channel had been formulated by combining two casted segments. The upper part and the lower part of the printed structure of the model were merged to form a microfluidic system to support a serpentine structure, for housing the printed construct. Media had been transfused through the serpentine structure from the inlet of the structure in the form of sinusoidal waveform and physiological waveform using a programmable syringe pump (DARWIN microfluidic Shenchen Syringe pump, Paris, France), as shown in Figure 4.

Uncertainty Analysis
Uncertainty analysis was performed to validate the model characteristics and to determine the shear stress sensitivity. The uncertainty equation of shear stress on the fabricated channel model was represented by Equation (5). Shear stress was generated when the printed serpentine model came into contact with the flowing HUVECS media.
Here, σ P , σ V and σ Y are the standard deviations of the pressure, velocity, and the strain rate standard deviation, respectively. From Equation (5), the uncertainty value of stress equivalent was evaluated as ±179.36501 N/m 2 .

Cellular Viability Assessment
The viability of the bioprinted vascular channel was evaluated over a period of 14 days by performing LIVE-DEAD assay. The channels were stained with the LIVE-DEAD cell imaging kit (Fisher Scientific, Waltham, MA, USA). Live cells were stained with Calcein-AM (2 µM), the fluorescence was measured at 494-517 nm, and dead cells were stained with ethidium homodimer-1 (4 µM), which was measured at 528-617 nm. The channels were incubated with LIVE-DEAD stain for 30 min and the samples were washed with PBS. Imaging was then performed using a Nikon C1 Confocal fluorescence microscope. Image J 1.53k open source software developed by Wayne Rasband and contributors National institute of Health, USA was used to evaluate the cell viability.

Sensitivity Analysis
The objective of the sensitivity analysis was to obtain the change of output with respect to change in input parameters. Input parameters may be a material property, variation in geometry and loading condition etc. For a proposed serpentine geometry, sensitivity analysis was performed for a no-slip wall condition of physiological flow corresponding to velocity model V 6 . The sensitivity of the serpentine model was calculated using Equation (6) where p is the calculated pressure, v is the maximum velocity (4.76 mm/min), ρ is the blood density 1066 kg/m 3 . Change in the pressure due to change in density can be expressed by the Equation (7) (modified version of Equation (6)) The analytical value of sensitivity was calculated by Equation (8) ∆p ∆ρ = 4.76 × 4.76 2 Sensitivity(analytical) = 11.32 (8b) The computed sensitivity was calculated from pressure obtained at blood density (1066 Kg/m 3 ) and 6% of blood density (997 kg/m 3 ) at a point P 1, respectively, as given in Equation The error of sensitivity was calculated between analytical sensitivity and computed sensitivity by below Equation (10) Error (sensitivity) = Sensitivity(computed) − Sensitivity(analytical) Sensitivity(analytical) Error (sensitivity) = 11.32 − 10.86 11.32 Error (sensitivity) = 4.06% At 6% of change in input parameter (density), it produces a minimum error of sensitivity [43]. From the above analysis it was concluded that an independent material property (density) of the flowing fluid influenced the dependent parameters (i.e., pressure) with a 6% change in density of fluid, as shown in Figure 5.

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 V1 and V3, there was no variation in the pressure parameter over time. Maximum pressure generation for velocity models V1, V2 and V3 were recorded as 2-5 × 10 5 Pa. Similarly, minimum pressure generation for different velocity models V1, V2, and V3 were obtained as 0.2 × 10 5 , −1 × 10 5 , and −1 × 10 5 Pa, respectively. In the velocity model V2, pressure generation on point P1 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 P1 and P3 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 (P1) and 20 mm/min in between the abdominal (P2) and rear regions (P3), 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 V2 at the end of the cycle in between the abdominal (P2) and rear regions (P3).

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 V1, V2, and V3 (Figure 7). A significant variation of pressure distribution was recorded for velocity model V1 compared to the two other velocity models. For velocity models V2 and V3, the local maximum pressure equal to 4 × 10 5 Pa was observed in the neck 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. 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 V1, V2, and V3 (Figure 7). A significant variation of pressure distribution was recorded for velocity model V1 compared to the two other velocity models. For velocity models V2 and V3, the local maximum pressure equal to 4 × 10 5 Pa was observed in the neck region (P1) 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 (P2). It was noticed that the value of axial velocity was similar for V1, V2, and V3 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 (P1) and periphery of the rear region (P3), 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 V1, V2, and V3 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 V4, V5, and V6 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 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.
Micromachines 2022, 13, x FOR PEER REVIEW 13 of 20 × 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 V4, 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 P1, 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 V4, V5, and V6 were obtained as 1 × 10 7 , 0.2 × 10 5 , and 0 × 10 5 , respectively. Similarly, 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 approx-imately 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. 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. 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.

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 P1, P2, and P3. 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

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.

Conclusions
A 3D bioprinting system using an extrusion-based technique was developed to fabricate vascular serpentine channels. The optimized perfusion rate of the media through the printed serpentine channels provided a biomimetic micro-niche for proliferation of seeded HUVECs. Dynamic cultured media at the extracellular surface of seeded cells imposed external forces on the membrane in the form of shear stress and its time derivatives on serpentine walls. Such a distributed pattern of extracellular forces induced faster and synchronous mechanotransduction to the focal points of adhesion of these cells. Hence, a variation in cellular proliferative activities was observed in response to the different perfusion rate of media. Therefore, it can be concluded that flow dynamics features, such as axial velocity, shear rate, stress, and pressure played a vital role in reproducing the physiological behavior and protruding adhesion properties of cells under in vitro conditions on the surface of the serpentine plane. Longitudinal stress generated by the flowing fluid also affects the HUVECS cell's functionality. It was also observed that the modeling of the serpentine with an elastic wall and consideration of physiological flow conditions at the inlet of the model brought maximal dynamicity in the axial velocity. On the other hand, identified local maxima and minima of flow parameters at designated regions of the serpentine were found to be a useful tool for effective seeding of cells for their proliferation over the surface of the serpentine. The minimum deviation of pressure, shear rate, and velocity were obtained for the sinusoidal flow with a slip wall condition. Result of Live-Dead assay showed the viability of HUVECs after 14 days of printing the channels, demonstrating that the obtained value of fluid dynamics parameters in the serpentine channel was helpful in maintaining the cell proliferation inside the bioreactor.