Numerical evaluation and experimental validation of fluid flow behavior within an organ-on-a-chip model

Computer Methods and Programs in Biomedicine 243 (2024) 107883 Available online 28 October 2023 0169-2607/© 2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Numerical evaluation and experimental validation of fluid flow behavior within an organ-on-a-chip model Violeta Carvalho a , b , c , d , e , * , Inˆ es M. Gonçalves b , f , g , Nelson Rodrigues c , Paulo Sousa d , e , Vˆ ania Pinto d , e , Graça Minas d , e , Hirokazu Kaji f , Su Ryon Shin a , Raquel O. Rodrigues d , e , Senhorinha F.C.F. Teixeira c , Rui A. Lima b , h , i a Division of Engineering in Medicine, Department of Medicine, Brigham and Women’s Hospital, Harvard Medical School, Cambridge, MA, 02139, USA b MEtRICs, Mechanical Engineering Department, University of Minho, Campus de Azur´ em, 4800-058 Guimar˜ aes, Portugal c ALGORITMI Center/LASI, University of Minho, Campus de Azur´ em, 4800-058 Guimar˜ aes, Portugal d Center for MicroElectromechanical Systems (CMEMS-UMinho), University of Minho, Campus de Azur´ em, 4800-058 Guimar˜ aes, Portugal e LABBELS—Associate Laboratory, Braga/Guimar˜ aes, Portugal f Institute of Biomaterials and Bioengineering (IBB), Tokyo Medical and Dental University (TMDU), 2-3-10 Kanda-Surugadai, Chiyoda, Tokyo 101-0062, Japan g Instituto Superior T´ ecnico, Universidade de Lisboa, Lisboa, Portugal h CEFT - Transport Phenomena Research Center, Faculty of Engineering, University of Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal i ALiCE - Associate Laboratory in Chemical Engineering, Faculty of Engineering, University of Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal ARTICLE INFO Keywords: Organ-on-a-chip CFD Computational simulations Microfluidics Biofluid mechanics Experimental validation ABSTRACT Background and Objective: By combining biomaterials, cell culture, and microfluidic technology, organ-on-a-chip (OoC) platforms have the ability to reproduce the physiological microenvironment of human organs. For this reason, these advanced microfluidic devices have been used to resemble various diseases and investigate novel treatments. In addition to the experimental assessment, numerical studies of biodevices have been performed aiming at their improvement and optimization. Despite considerable progress in numerical modeling of biodevices, the validation of these computational models through comparison with experimental assays remains a significant gap in the current literature. This step is critical to ensure the accuracy and reliability of numerical models, and consequently enhance confidence in their predictive results. The aim of the present work is to develop a numerical model capable of reproducing the fluid flow behavior within an OoC, for future investigations, encompassing the geometry optimization. Methods: In this study, the validation of a numerical model for an OoC microfluidic device was undertaken. This comprised both quantitative and qualitative assessments of trace microparticles flowing through a physical OoC model. High-speed microscopy images of the flow, using a blood analog fluid, were analyzed and compared with the numerical simulations run using the Ansys Fluent software. For a qualitative analysis, the particles’ paths through the inlet and bifurcations were observed whereas, for a quantitative analysis, the particle velocities were measured. Furthermore, oxygen transport was simulated and evaluated for different Reynolds numbers. Results: In both qualitative and quantitative analyses, the results predicted by the numerical model and the ones outputted by the experimental model were in good agreement. These findings underscore the capability and potential of the developed numerical model. The examination of oxygen transport at various vertical positions within the organoid has revealed that for lower positions, oxygen transport predominantly occurs through diffusion, leading to a symmetric distribution of oxygen. Contrastingly, the convection phenomenon becomes more evident in the upper region of the organoid. Conclusions: The successful validation of the numerical model against experimental data shows its accuracy and reliability in simulating the fluid flow within the OoC, which consequently can expedite the OoC design process by reducing the need for prototypes’ fabrication and costly laboratory experiments. * Corresponding author. E-mail address: [email protected] (V. Carvalho). Contents lists available at ScienceDirect Computer Methods and Programs in Biomedicine journal homepage: www.elsevier.com/locate/cmpb https://doi.org/10.1016/j.cmpb.2023.107883 Received 8 August 2023; Received in revised form 22 October 2023; Accepted 23 October 2023 Computer Methods and Programs in Biomedicine 243 (2024) 107883 2 1. Introduction Drug development and consequent approval for clinical use is a timeconsuming and costly process, and often ineffective [1–3]. This is usually due to the lack of translatability of current preclinical models to evaluate the effectiveness and toxicity of new drug candidates [4–6]. In response to this issue, three-dimensional (3D) in vitro models such as organ-on-a-chip (OoC) have been progressively used in biomedical research and have contributed to ground-breaking progress in life sciences [7–11]. OoC technology allows achieving an in-depth understanding of the disease mechanisms and testing novel therapies by recreating the dynamic microenvironment of human tissues and, consequently, increasing their efficiency [12–14]. In addition, the 3Rs principle which looks for Replacing, Reducing, and Refining the use of animal models in research, can be followed [15–17]. Computational tools offer several advantages over experimental methodologies since they provide an overview of flow physics with good precision and accuracy in a rapid and cost-effective way [18–23]. The use of these tools has been increasing in both industry and research of wide engineering fields, namely in microfluidics. This field of research has rapidly evolved and the scientific community has combined computational tools alongside theoretical and experimental methods [24–29]. However, an important task to perform before having confidence in the numerical model relies on its validation through empirical tests, both qualitatively and quantitatively [30–35]. This is usually achieved through measurements of velocity [36–39], pH [40], temperature [41–43], droplet formation [44,45] and length [46], shear stress [47], colorimetric techniques [48], particles width and path [49–51], mixing quality [52–55], potential and power density [56,57], impedance [58], trapping process [59], among others performed by the researchers themselves, or, in other cases, the numerical model is validated with data published by other researchers [53,60–68]. For instance, since reactive oxygen species, such as hydroperoxide (H 2 O 2 ), are known to be a major factor in cancer development, Lo et al. [69] evaluated the concentration of H 2 O 2 , both numerically and experimentally, and the results showed to be in good agreement. On the other hand, Komen et al. [70] developed a microfluidic device where cells were exposed to an in-vivo-like drug concentration profile and this environment was also simulated through computational tools. For validation, the authors conducted fluorescence assays and the results were similar to the numerical ones. Mosavati and the research team [71] evaluated, numerically, the convection-diffusion phenomena of glucose through a porous membrane and validated it by comparing, experimentally, the measured glucose concentration at different positions in the microfluidic device. Numerical models can also be validated through velocity measurements. For instance, Gracka and co-workers [39] developed a numerical model of a hyperbolic microchannel able to simulate the cell-free layer effect considering a multiphase approach. The results were validated through videos recorded during the experiments. A manual tracking method was used to follow the red blood cells individually and then their velocity was estimated. Similarly, Zhang and co-workers [36] simulated a multiphase flow within a cells-on-a-chip model and validated it by tracking tracing particles using the microscope. A good qualitative and quantitative agreement was obtained between the experimental and simulation results for particle velocity in the microfluidic chip. On the other hand, Tanaka et al. [37] used image velocimetry (PIV) technique to validate the computational fluid dynamics simulation of an autonomous hybrid pump driven by self-organized microtissues of cardiomyocytes. The same technique of velocity measurements was used by Hynes et al. [38]. The authors investigated computationally the fluid flow behavior considering both acellular and endothelialized branching vessels. The computational results showed to be in good accordance with the experimental measurements both qualitative and quantitatively. Other authors consider the numerical model validated since the experimental outcomes met the expectations, but they do not perform a graphical comparison of both models [72,73]. While numerical simulations for OoC research have been progressively increased, validation of OoC numerical studies remains scarce in the literature. Besides, to the best of our knowledge, it is the first time that an OoC numerical simulation has been validated by evaluating its fluid flow and respective velocity magnitudes by means of a particle-tracking velocimetry method. Beyond validation, the present work offers new insights through in-depth fluid flow analysis. It should be also noted that current OoC development often relies on trial and error when setting the flow rates and other cell culture parameters, hindering its progress. The present study, on the other hand, showcases numerical simulations’ potential to expedite optimized OoC systems, validated experimentally. By emphasizing synergy between modeling and experimentation, its systematic use is promoted, advancing OoC technology and reducing prototypes and lab costs. 2. Methods In this section, the device fabrication and experimental setup are first presented. Then, the respective numerical model as well as the simulation details namely the mesh, governing equations, boundary conditions, and solver, are detailed. 2.1. Microfluidic device fabrication The OoC microfluidic device was fabricated in polydimethylsiloxane (PDMS) by soft lithography technique, using a mold of polymethyl methacrylate (PMMA) manufactured by laser cutting. PDMS was chosen for this study due to its unique properties and suitability for the research objectives, namely biocompatibility, gas permeability, ability to culture cellular models, transparency, and easy/cheap fabrication. The good Fig. 1. (a) Top view of the PDMS OoC microfluidic device used to perform flow visualizations and (b) three-dimensional geometry of the microfluidic device and its main parts, the inlet, outlet, and culture compartments. V. Carvalho et al. Computer Methods and Programs in Biomedicine 243 (2024) 107883 3 optical transparency of PDMS is particularly important as it allows to perform the experimental flow visualizations and in this way, it is possible to compare with the numerical simulations. [74]. After mold fabrication, 25 g of a PDMS solution, composed of 1/10 %wt of base/- curing agent (Sylgard® 184 Silicone Elastomer Kit) was poured over the PMMA mold and cured at 80 ◦C for 45 min. After the PDMS curing, the structure was removed from the mold and holes were then punched to provide access to the inlets and outlets, placed in desired locations. Finally, the channels were closed with a glass slide by oxygen plasm for 30 s at 30 W (Plasma Systems ZEPTO from Diener electronic). The PDMS microfluidic device, shown in Fig. 1, was developed by our group for cellular studies of a set of organoids (i.e., to contain four biomodels inserted in parallel). Nevertheless, in the present work, the microfluidic chip was used for flow validation purposes, and thus, without cells/organoids. 2.2. Blood analog and experimental setup for flow visualizations For the experimental studies, the fluid velocity was evaluated by the tracking of particles from a blood analog fluid flowing through the OoC microfluidic device. The blood analog fluid was prepared by mixing the surfactant Brij L4 (Merck KGaA, Darmstadt, Germany) in distilled water at a concentration of 1 wt.% with a final density of 996 kg/m3 and a viscosity of approximately 0.0012 Pa ⋅ s [75]. The solution was then mechanically stirred until became homogeneous. Lastly, the solution was pushed two times through a precolumn with a filter with pores of 20 μ m in size. This step allows more uniformity of the formed micelles, as well as guarantees that their dimensions are similar to, or smaller than, the pore size of the filter. The OoC microfluidic device was placed in an inverted microscope (IX71; Olympus Corporation, Tokyo, Japan) connected to a high-speed camera (Fastcam SA3, Photron, Motion Engineering Company, Westfield, IN, USA). The blood analog fluid was loaded into a syringe placed on a syringe pump (NEMESYS, CETONI, Germany) and connected to the OoC microfluidic device, as schematized in Fig. 2a). Firstly, the flow rate of the fluid was set on the pump to 600 μ L/min to measure velocities at the center of one compartment of the device and qualitatively validate the flow streamlines at the inlet, and bifurcations. For additional qualitative validation, higher flow rates were tested to evaluate the variations in flow behavior including 1000, 1500, and 2000 μ L/min, as well as a lower flow rate of 336 μ L/min. Regardless of the flow rate, all videos of the flow were acquired at a frame rate of 3000 frames/s. The particle velocity measurements were conducted by using the open-source image analysis software ImageJ, where a manual tracking plugin, MTrackJ, was used to track the particles individually as shown in Fig. 2b). By using this plugin, the particles’ velocity was calculated by determining the position of the selected in-focus particles in several consecutive frames of the recorded movie. For each plane, only in-focus particles were selected, while out-of-focus particles were not considered. Regarding the velocity calculation, this was based on the distances that a particle was moved for a knowing time according to the selected frame rate of the high-speed camera. By calculating the particle flow velocities at different planes, it was possible to determine the velocity profile at the culture compartments. It should be noted that the geometry is assumed Fig. 2. (a) Experimental setup constituted by an inverted microscope connected to a high-speed camera and a syringe pump used to visualize the fluid flow, (b) Scheme of the region of the particle tracking process, at the central region of the culture compartments and indication of the position of bifurcation 1 (B1) and bifurcation 2 (B2). V. Carvalho et al. Computer Methods and Programs in Biomedicine 243 (2024) 107883 4 as symmetric, and thus both the velocity measurements obtained experimentally and the velocities exported from Ansys were analyzed only in one of the four compartments. 2.3. Microfluidic device geometry and mesh The geometry and dimensions of the developed microfluidic device are represented in Fig. 3. Note that the height of the OoC microfluidic device is 5 mm. The first step for the numerical studies consisted of developing three different hexahedral meshes with 156,372; 325,068, and 683,370 elements. Aiming to verify the precision of the simulations, the grid convergence method (GCI) was used to analyze the mesh quality and the discretization error through the evaluation of the velocity at Z =1 mm [76]. The GCI is calculated through Eq. (1): GCI21 fine =1.25e21 a rp 21 −1(1) e21 a= ∅1−∅2 ∅1 (2) p=1 ln(r21) ln ε 32 ε 21 (3) Being r 21 the grid refinement factor, e21 a the approximate relative error (Eq. (2)) and p the apparent order (Eq. (3)), with ∅ 1 and ∅ 2 the velocity computed with the finest and medium mesh, respectively, and ε 32 =∅ 3 −∅ 2 . It should be noted that indexes 21 and 32 represent the comparison between the finest and medium mesh, and between the coarse and medium mesh, respectively. The GCI decreased substantially from 0.32 % to 0.04 % (Fig. 4a)), i. e., with further refinement the computed velocities would not change considerably, and thus, the medium mesh, with 325,068 elements was used in numerical simulations (Fig. 4b)). 2.4. Governing equations The fluid flow through the microfluidic channels was simulated by solving the Navier–Stokes equations for mass and momentum continuity Eqs. (4) and ((5)): ∇v →=0(4) ∇⋅( ρ v →v →) = − ∇p+ μ ∇2v →+ ρ g →(5) where v →is the velocity vector, p is the static pressure, ρ is the fluid density, μ is the dynamic viscosity, and g →the gravity vector. Given the high flow rate used to run the experimental assays, the transitional shear stress transport (SST) kω turbulence model was also tested for validation purposes in case of turbulent flow occurs experimentally. This is one of the Reynolds-averaged Navier–Stokes (RANS) models available in Ansys Fluent software and it consists of a set of timeaveraged equations with turbulent features (turbulent kinetic energy, k, and the specific dissipation rate, ω Eqs. (6) and ((7))) [77,78]. ∂ ∂ xi ( ρ kui) = ∂ ∂ xj(Γk ∂ k ∂ xj)+Gk−Yk+Sk(6) ∂ ∂ xi ( ρω ui) = ∂ ∂ xj(Γ ω ∂ ω ∂ xj)+G ω −Y ω +D ω +S ω (7) where, G k represents the generation of k and G ω represents the generation of ω ; Γ ω and Γ k represent the effective diffusivity of k and ω , respectively; Y k and Y ω correspond to the dissipation of k and ω due to turbulence; D ω represents the cross-diffusion term and S k and S ω are userdefined source terms. After validating the numerical model, the presence of hydrogel scaffolds in the culture chambers was considered. These are threedimensional culture models used to mimic organ functions. In this case, these were considered solid bodies, an approximation previously considered by other authors [79], and the oxygen transport was simulated considering an arbitrary scalar (ϕ). For steady-state simulations, ANSYS Fluent will solve Eq. (8): ∂ ∂ xi( ρ uiϕ−Γ ∂ ϕ ∂ xi)=S∅(8) where Γ and S ∅ are the diffusion coefficient and source term. 2.5. Boundary conditions The walls of the microfluidic chip were considered rigid with a noslip condition and a zero-gauge pressure condition was set at the outlet of the device. At the inlet, a constant flow velocity was imposed which matched the flow rates used on the flow visualization tests. Regarding the fluid properties, for validation purposes, the numerical simulations were run considering the properties of the blood analog Fig. 3. Geometry and dimensions (in millimeters) of the used microfluidic device. V. Carvalho et al. Computer Methods and Programs in Biomedicine 243 (2024) 107883 5 (a) (b) 0.0 0.1 0.2 0.3 0.4 GCI32 GCI (%) GCI21 Fig. 4. (a) Mean percentage of the grid convergence index - GCI obtained from the computed velocities at Z =1 mm between medium and fine mesh (GCI21) and coarse and medium mesh (GCI32); (b) Hexahedral mesh used in the numerical simulations with 325,068 elements. V. Carvalho et al. Computer Methods and Programs in Biomedicine 243 (2024) 107883 6 fluid, BrijL4, used experimentally, assuming a single-phase approach. On the other hand, for the remaining computational simulations, the presence of cells within hydrogel scaffolds modeled as solid bodies was considered, as well as the properties of the culture media used in cellular assays ( ρ =1007 kg/m3 and µ=0.958 ×10–3 kg/(m.s)). In addition, oxygen transport was simulated taking into consideration the fluid closed circuit, the oxygen consumption by cells, and different oxygen diffusion in both the hydrogel scaffold and the culture media [22,23,80, 81]. Herein, the hydrogel scaffold was considered to have a density of 150,000 cells. 2.6. Solver ANSYS Fluent software was used to solve the Navier-Stokes equations through the Finite Volume Method approach. In this method, the domain is divided into control volumes, and then the equations are solved. The values of the various flow variables are calculated at the centroid of each control volume. In this study, the pressure-based coupled algorithm was used since it is considered adequate for steadystate simulations and the convergence was attained when the residuals were inferior to 10 −6 for all simulations. Flow Direction Numerical Experimental (a) (b) (c) B1 B2 Fig. 5. Qualitative validation through velocity streamlines from numerical simulations (left) and Z-project image from experimental assays (right) at (a) inlet, (b) B1, and (c) B2) for the flow rate of 600 µL/min. Flow direction Fig. 6. Three-dimensional fluid flow along the OoC microfluidic device with inlet region zoomed to verify the vortices’ creation. V. Carvalho et al. Computer Methods and Programs in Biomedicine 243 (2024) 107883 7 3. Results In this section, the validation of the numerical model is presented, firstly qualitatively by comparing the fluid flow behavior in terms of streamlines, and then quantitatively by comparing the particles’ velocity measured experimentally with the ones predicted by the numerical model. In addition, a detailed analysis of the fluid flow is presented. 3.1. Qualitative validation As the first step for validation, a qualitative analysis of the flow was conducted. For this purpose, the Z-project tool of ImageJ software was used and the obtained images from the experimental assays were compared with the velocity streamlines observed in Ansys Fluent software at Z =2 mm. The comparison made for the flow rate of 600 µL/min is depicted in Fig. 5 for the inlet (Fig. 5a, Video S1), bifurcation 1 (Fig. 5b, Video S2) and bifurcation 2 (Fig. 5c, Video S3). Fig. 5a) presents the streamlines at the inlet region. With the selected flow rate (600 µL/min), vortices are created as can be observed in both numerical and experimental images. Fig. 5b) and c) present the comparison of the flow behavior in the first two bifurcations of the OoC microfluidic device, whereas B1 has a higher angle than the second one (B2), as shown in Fig. 2b). Despite the creation of recirculation zones in the inlet, the alignment of the particles with the flow can be clearly observed at both bifurcations. In all cases, the numerical and experimental results are in close agreement. Since the previous results were obtained from an XY plane, in Fig. 6, the three-dimensional fluid flow in the three regions is represented with normalized velocity vectors. Here, the recirculation at the inlet is more perceptible. Also, it can be seen that immediately after the inlet region, the flow is aligned. In addition to the previous experiments, different flow rates were tested for qualitative validation as shown in Fig. 7. The flow differences were noticed only at the inlet. For the lowest flow rate (336 µL/min), no vortices are observed both experimentally and numerically. Nevertheless, by increasing the flow rate, the creation of vortices increased at the device entrance. Despite this, at both bifurcations, the flow is fully developed and laminar. These data are presented in supplementary materials. 1500 μL/min 2000 μL/min Flow Direction Numerical Experimental 336 μL/min 1000 μL/min Fig. 7. Qualitative validation of velocity streamlines from the numerical simulations and through Z-project images from the experimental assays, at the inlet for different flow rates (336, 1000, 1500, and 2000 μ L/min). 0246 0.00 0.05 0.10 0.15 0.20 Experimental Numerical (Laminar) Numerical (k-Z SST model) ZPosition (mm) Velocity (mm/s) Fig. 8. Velocity profile obtained through numerical simulations considering both laminar and turbulent (kω SST model) approaches, and experimental data obtained from particle tracking. Error bars represent the standard deviation. V. Carvalho et al. Computer Methods and Programs in Biomedicine 243 (2024) 107883 8 3.2. Quantitative validation To complement the previous results, a quantitative validation was carried out by measuring, experimentally, the velocities in the culture chambers, since they are the region of interest, where cell-derived organoids would be placed. For these experiments, only the 600 μ L/ min flow rate was used, since this was the flow rate that allowed performing the particle tracking with more accuracy. To this end, particles’ velocities were measured experimentally from the upper to the bottom of the channel at different heights along the Z-axis of the circular compartments (Fig. 2b)). The comparison was made with the velocities computed along the line schematized in Fig. 8. As previously mentioned, both laminar and turbulent (kω SST model) approaches were evaluated in quantitative terms to assess if the predicted velocity magnitude would vary. The respective velocity profiles, as well as the experimental data, are depicted in Fig. 8. A good agreement between the numerical and experimental results is observed, mainly near the walls where the velocities are lower. In the central region of the parabolic profile, the measured experimental velocity is slightly higher than the numerically predicted. Since in this region, particles have higher velocities, motion blur occurs and the measurements become less accurate [78,82]. On the contrary, near the Fig. 9. Computed 3D velocity profiles in the microfluidic device through CFD results at different YZ planes along the domain. Fig. 10. Velocity streamlines and wall shear stress around and at the organoids for a) Re =0.26, (b) Re =0.47, and (c) Re =1.56. V. Carvalho et al. Computer Methods and Programs in Biomedicine 243 (2024) 107883 9 walls, the velocities are lower, and the tracking results are more accurate. Furthermore, one can see that the velocity profiles for laminar and turbulent models overlapped. Despite this being visible, SPSS software was used to evaluate the existence of statistically significant differences (t-test for dependent variables, p >0.05). Since no significant differences were detected, it can be concluded that despite vortices being created at the inlet of the device, the flow develops as laminar along the microfluidic channels, therefore, the laminar model is adequate to accurately predict the velocity magnitudes. 3.3. 3D velocity profiles After validating the numerical model, the CFD three-dimensional velocity profiles were computed considering the same flow rate as previously (600 μ L/min) and evaluated along the device focusing not only on the culture chambers, but also on the different bifurcations sites B1 and B2 (Fig. 9). Looking at the previous results, one can see that in plane A (B1), the flow is about to separate. This phenomenon can be better examined further ahead in B2. Here, three different planes were evaluated, and the evolution of the fluid flow can be clearly observed. In plane B, the flow is fully developed with a parabolic profile, and as it approximates the bifurcation, the velocity profile changes its shape (plane C) until fully separated (plane D), creating two peaks in the plot due to the symmetry of the geometry. On the other hand, before the culture chamber (plane E), the velocity profile is again parabolic. Although at the culture Fig. 11. Oxygen profiles at XY planes for (a) Z =0 mm, (b) Z =1 mm, and (c) Z =2 mm, which corresponds to the bottom, middle, and top positions of the organoid for different Reynolds numbers Re=0.26, Re=0.47, and Re=1.56. V. Carvalho et al.