Accurate flow in augmented networks (AFAN): an approach to generating three-dimensional biomimetic microfluidic networks with controlled flow

In vivo, microvasculature provides oxygen, nutrients, and soluble factors necessary for cell survival and function. The highly tortuous, densely-packed, and interconnected three-dimensional (3D) architecture of microvasculature ensures that cells receive these crucial components. The ability to duplicate microvascular architecture in tissue-engineered models could provide a means to generate largevolume constructs as well as advanced microphysiological systems. Similarly, the ability to induce realistic flow in engineered microvasculature is crucial to recapitulating in vivo-like flow and transport. Advanced biofabrication techniques are capable of generating 3D, biomimetic microfluidic networks in hydrogels, however, these models can exhibit systemic aberrations in flow due to incorrect boundary conditions. To overcome this problem, we developed an automated method for generating synthetic augmented channels that induce the desired flow properties within three-dimensional microfluidic networks. These augmented inlets and outlets enforce the appropriate boundary conditions for achieving specified flow properties and create a three-dimensional output useful for image-guided fabrication techniques to create biomimetic microvascular networks.


Introduction
The ability to generate vascularized tissue constructs has been a major challenge in the eld of tissue engineering for decades. [1][2][3] To overcome this limitation, small volume constructs have been implemented that rely on diffusionmediated transport to deliver oxygen, nutrients, and soluble factors, and to remove waste products. The ability to generate hydrogels laden with microvasculature that recapitulates the dense and tortuous architecture of in vivo vascular networks could aid in fabricating larger volume tissue constructs 4 and advanced cell culture platforms to model physiological and pathological processes.
Two general approaches to generate vascularized constructs have been developed: self-assembly of networks by vascular cells and formation of pre-dened, engineered vasculature followed by vascular cell seeding. Self-assembly provides a straightforward approach to create microvasculature in both natural [5][6][7][8][9] and synthetic 10,11 hydrogels. When performed in hydrogels housed in a microuidic device, self-assembled networks can anastomose with larger microuidic channels allowing for uid ow through the networks. 5,8,12 This approach has been implemented to investigate leukocyte adhesion to vessel walls, 8 cancer cell extravasation, 13 convective transport, 10 soluble factor signaling during vasculogenesis, 14 and endothelial cell response to uid ow. 8 While self-assembly is a powerful approach to generate microvascular networks, there is no control over the nal architecture, therefore making it impossible to use the same network across multiple experiments or to easily couple experimental data with computational uid dynamics (CFD) models. Furthermore, self-assembly has not been implemented to generate larger diameter arterioles and arteries needed to reproduce the hierarchical structure of in vivo vasculature.
To overcome these limitations, various biofabrication techniques have been developed that allow pre-denition of network architecture prior to seeding with vascular cells. Many of these techniques are amenable to generating twodimensional (2D), planar networks embedded in hydrogels which have been implemented for drug screening, 15,16 blood ow modeling, [17][18][19] disease modeling, 20,21 to investigate owmediated signaling between tumor and endothelial cells, 22 cancer cell extravasation, 23 self-healing of vasculature post inammation, 24 and thrombotic response of vessels due to inammation. 20 While many existing microfabrication approaches allow for direct control over 2D, planar networks, they do not duplicate the complex 3D structure of in vivo networks. 25 A few techniques to generate 3D, non-planar networks have been developed, including 3D printing of sacri-cial carbohydrate glass, 26 modular assembly, 27 and 3D bioprinting. 28 These approaches have been implemented to increase cell viability in larger volume constructs, 26 to create a perfusable microuidic hydrogel, 27 and to spatially organize both cells and vasculature in a tissue construct. 28 Although these approaches allow repeated fabrication of 3D vascular networks with well-dened geometry, they are incapable of fabricating dense and tortuous small-diameter structures needed to recapitulate in vivo microvascular architecture. Direct-write assembly 29 and omnidirectional printing, 30 are capable of generating a hierarchical vascular system with a wide range of diameters, 10-530 mm for direct-write assembly and 18-600 mm for omnidirectional printing, but are not yet amenable to generating dense, tortuous, and highly interconnected microvascular networks. To overcome this problem, we recently developed an image-guided, laser-induced hydrogel degradation (LIHD) technique that utilizes either computer aided design (CAD) synthetic networks or 3D image stacks of in vivo vasculature as digital templates to fabricate 3D, biomimetic, hydrogel-embedded microuidic networks whose architecture closely matches that of in vivo microvasculature. 31,32 Laser-induced degradation is amenable to both synthetic and natural hydrogels [31][32][33][34][35][36][37][38][39] and has been utilized to generate micro-uidic networks in cell-laden constructs. 40 One limitation of using 3D image stacks of in vivo vasculature as a digital template is the presence of dead-end vascular structures that impede ow (Fig. 1c). Similar to modeling synthetic microvascular networks (SMNs), real microvascular networks (RMNs) must be cropped to t within designed volumes. Since RMNs exhibit a high degree of connectivity, the selection of any region results in an incomplete network with multiple terminations, making it difficult to prescribe desired boundary conditions. Two generalized methods have been proposed to overcome this problem: incorporating multiple inlets and outlets to the network (Fig. 1b) 41 or completing the network with additional connections to force the terminating vessels to converge to a single bounding channel (Fig. 1a). 40 While the rst approach is capable of generating physiologically realistic ow through the network, 42 it is limited to 2D, planar networks. This approach also requires the implementation of multiple syringe pumps and connectors which is experimentally cumbersome. The second method (Fig. 1a) offers a simple solution to this problem by adding an additional bounding channel that connects all of the network dead-ends to a single inlet and outlet, but has only been implemented for 2D, planar networks and does not provide control over which dead-ends are inlets or outlets, thereby limiting control over ow properties. To overcome these problems, we developed a network design approach that constructs synthetic connections between boundary nodes and two feeders (Fig. 1d), an inlet and outlet, allowing for well-controlled boundary pressure, at previous dead-ends based on established models or in vivo measurements. Accurate ow simulations are required to calculate these boundary connections. Rigorous numerical methods, such as the nite element method 43 and nite difference method 44 are frequently used in CFD. However, determining the parameters for the necessary augmented elements would require a timeconsuming iterative optimization step that is impractical on most workstations. Accordingly, we performed this optimization using the linear Hagen-Poiseuille (H-P) method, which allows for a rapid extrapolation of pressure-driven ow through fabricated networks. We implemented these features in an interactive soware package called Accurate Flow in Augmented Networks (AFAN) using C++ and CUDA for fast evaluation and real-time visualization.
We demonstrate that the resulting augmented network design can be fabricated using LIHD. We also demonstrated that our augmented microvascular model has a high level of ow predictive capability based on comprehensive CFD models and microuidic experiments. This approach lays the foundation for implementing 3D image stacks as digital templates for fabrication of biomimetic vascular networks embedded in hydrogels whose architecture and ow properties accurately mimic in vivo vasculature.

Two-dimensional Wheatstone bridge network
We designed and fabricated a microuidic network based on the Wheatstone bridge using the AFAN user interface. Based on user input describing network architecture and ow properties, AFAN generates a binary mask used as the basis for fabrication.
The fabricated network was constructed based on our previously published protocols. 32 Plasma-bonded polydimethylsiloxane (PDMS) and glass microuidic devices were fabricated to provide connections for a syringe pump. For hydrogel incorporation within a device, devices were functionalized by washing with 2,2-dimethoxy-2-phenylacetophenone, and 3-(trimethoxysilyl)propyl methacrylate, to allow hydrogel bonding to the PDMS and glass. 10 A pre-polymer solution of 5% weight per volume 3.5 kDa poly(ethylene glycol) diacrylate (PEGDA), 3.7 mM Alexa Fluor 633-labeled acryl-poly(ethylene glycol), and 10.2 mM lithium phenyl-2,4,6-trimethylbenzoylphosphinate (LAP) in HEPES-buffered saline (HBS) (pH 8.3, 10 mM HEPES, 100 mM NaCl) was photopolymerized inside the microuidic device via exposure to UV light at 6 mW cm À2 for 30 seconds. A photo-mask placed in the light path was used to dene the geometry of the hydrogel within the device. LIHD was used to create a microuidic network in the PEGDA hydrogels as previously described. [31][32][33] A series of virtual masks 45 dening ROIs in x, y, and z, were generated to guide the position of a 140 fs pulsed Ti:S laser operating at 790 nm at 37.7 nJ mm À2 focused through a 20Â (NA ¼ 1.0) water immersion objective for selective hydrogel degradation.
The fabricated network was lled with 2000 kDa FITClabeled dextran at 1 mg mL À1 in HBS and imaged via structured illumination. Particle image velocimetry was used to quantify the average velocity in each segment of the micro-uidic network. Using a syringe pump, 3 mm diameter polystyrene spheres, at 8.4 Â 10 6 particles per mL in HBS, were owed through the network at 25 mL h À1 . Images of the particles were acquired using a 2 ms image acquisition time over 3 min intervals. An average of 200 particles per 3 min interval were analyzed per segment in triplicate. The center-to-center distance traveled by each particle was measured using ImageJ and the particle velocity was calculated by dividing the distance traveled by the 2 ms acquisition time. Particles that spanned segments, or that overlapped, were excluded and particle velocities were averaged for each segment. Constant uid ow was maintained and veried throughout the image collection period. The particle streaks collected are colored to help distinguish particles from each other (Fig. 6b).

Three-dimensional microvascular network
The overview for our approach is shown in Fig. 2. We used a whole mouse brain microvascular data set (Fig. 7a) collected using knife-edge scanning microscopy (KESM) 46,47 available through the KESM Brain Atlas (https://www.kesm.cs.tamu.edu). 48 A 652 Â 652 Â 100 pixel (120 Â 120 Â 100 mm) region of interest (ROI) was identied and extracted from the whole-brain data set. Microvessel centerlines and connectivity were segmented using a predictor-corrector algorithm, 49 while the surface model and radii were extracted manually by setting a threshold to separate the microvascular structure from the background. This data was combined to generate a graph-based model used as input to the AFAN soware.
The AFAN user interface was used to specify boundary pressures producing ow velocities based on in vivo measurements. 50,51 The AFAN H-P method was used to extrapolate ow properties throughout the ROI. The segmented network was then augmented with connections to enforce the desired boundary conditions. The augmented network was constructed such that an input volumetric ow rate Q i produced the desired ow characteristics in the ROI (Fig. 1d).
A comparative CFD model was generated from an AFANproduced binary mask and imported into OnShape (https:// www.onshape.com), an online CAD package. The geometric mesh was constructed using a built-in meshing function in SimScale (https://www.simscale.com), which creates polyhedral meshes for uid-based models. SimScale was used to integrate an incompressible steady laminar uid ow into the simulation with the simpleFoam solver in the OpenFOAM toolbox. 52 We specied the appropriate materials and boundary conditions (inlet velocity, outlet pressure, and no slip walls), and conducted simulations on the cloud. The results were visualized and analyzed using ParaView, 53 which offers a host of postprocessing operations for data analysis. The augmented, 3D biomimetic network was fabricated in PEGDA hydrogels and visualized with uorescent dextran in the same manner as the 2D network.

Network characterization
Characterizing the hydraulic resistance of each synthetic connection relies on accurate boundary measurements of either boundary pressures or ow rates. The proposed interpretation of pressure-driven ow through circular microchannels uses the H-P equation based on boundary pressures: which estimates the volumetric ow rate Q as a function of the pressure drop DP, dynamic viscosity m, channel length L, and channel radius r. If necessary, the corresponding wall shear stress s can be calculated from the computed average ow velocity v based on eqn (2).
The H-P equation was initially derived for channels that are innitely long with no variation in geometry, but microvessels oen have nite lengths and changing cross-sections. Since the H-P method builds on several assumptions, including Newtonian, incompressible, and laminar ow properties, and boundary conditions such as the uniform pressure gradient condition, we provide the following justication to demonstrate its viability in our proposed models: for ow regimes, in vivo blood ow exhibits a Newtonian pattern when the shear rate is greater than or equal to 100 s À1 , 54 and small microvessels such as capillaries and arterioles satisfy this condition. 55 Furthermore, uids in capillaries are well approximated as incompressible mediums as the Mach number Ma drops below 0.2. 56 Finally, the Reynolds number Re in capillaries is usually much smaller than 2000, and oen less than 1, therefore such ow is completely laminar. As for boundary conditions, the H-P equation can be reasonably applied to channels with nite lengths if the uid has a local and fully developed laminar ow prole. 57 While a modi-ed H-P equation was developed to correct for the tortuous and fractal properties of capillaries, 58 a capillary ber can also be approximated as a straight cylinder of the same length with no variation in diameter since inertial effects are negligible and the cross-sectional velocity is stable for curved paths in low Reynolds situations. 59 In conclusion, the H-P method can be applied to extrapolate accurate pressure-driven ow in our proposed models.
Nevertheless, eqn (1) only quanties the average velocity through one channel, requiring an integral method that links all channels together to solve for the entire network. A common technique is to apply an electric-hydraulic analogy using the following strategy: 57 (1) Transform the hydraulic network into an equivalent electric circuit and complete the circuit with voltage sources (resembling pressure sources) based on an educated guess or previous measurements.
(3) Write an Ohm's function for each ber and a Kirchhoff's current function 57 for each branching point.
(4) Organize these equations into a linear system and solve with standard matrix factorization.
The high-dimensional matrix factorization step can take hours using CPU-based implementation for large networks. In this work, this step was performed using the CUBLAS library in combination with an nVidia Geforce 970 GTX graphics card to provide real-time performance, and interactive visualization and renement of the network and its augmented components.

Network augmentation
Since modeling entire microvascular networks is experimentally impractical, current studies rely on selected ROIs that are culled from a larger network. The peripheral network of an ROI plays an important role in regulating its boundary pressure. Simply building a microuidic model consisting of an extracted ROI will not produce predictable ow as most boundary nodes are either open to atmosphere or abruptly terminated, and multiple source supplies can be attached to provide the necessary boundary conditions. However, this method is impractical for 3D networks that have vessels terminating from all directions of the ROI volumes. We addressed this by developing a network augmentation technique that connects all the inlet nodes to an inlet feeder and all the outlet nodes to an outlet feeder using synthetic connections. Not only do these connections compensate for the absence of the peripheral network in controlling boundary conditions, but they also reduce the number of sources connected to a microuidic device and therefore can be easily integrated into most microuidic devices.
The augmented network is composed of three parts: an ROI network, augmented connections, and two main feeders. The synthetic connections play the role of a ow splitter to simultaneously achieve the appropriate boundary pressure P i at each end while keeping the ROI network unchanged (Fig. 1d). Since augmenting outlets is identical to inlets, we will focus our discussion on inlets in the following sections.
According to eqn (1), building a synthetic connection with the appropriate hydraulic resistance requires a trade-off between channel length and radius. For planar networks, we xed channel length and solved for radius as an easy way to avoid intersection (Fig. 3a). There are several graphical user interface libraries, such as the OpenGL library, that simplify 2D spatial arrangements. To build xed-length connections, we introduced a pool-like feeder (a cylinder) with uniform pressure outputs around the surface. To obtain the source pressure, the inlet corresponding to the lowest source pressure was constrained to 5 mm. We expanded other connections, decreasing their resistance to match the calculated source pressure.
For non-planar networks, we introduced a bus-like feeder (a cuboid) with uniform pressure outputs at the top and bottom faces (Fig. 3b), which connes connections to a single plane (e.g. x-y). This helps determine optimal arrangements while avoiding overlaps. To obtain the source pressure, we set all radii to 5 mm and xed the length of the inlet requiring the highest source pressure. Other connections were extended to increase their resistance to satisfy the established boundary conditions, or calculated source pressures. Unlike channel expansion, channel extension must account for channel distribution and space utilization. While multiple methods can be used to increase the channel length, one can simply replace the original connections with longer straight connections. Unfortunately, this approach sacrices space efficiency for simplicity and may make the network difficult to fabricate. More complex paths, such as Hilbert curves, 60 provide an efficient use of space, but are incompatible with our fabrication method as our hydrogels may not be able to mechanically support such a densely-packed structure. We address this by adopting a space-lling method using square curves (Fig. 4), which provides a compromise between simplicity, space efficiency, and fabrication integrity. Each of the square wave-like connections was constructed inside of a user-dened axis-aligned bounding box (AABB), which allowed for rapid testing of intersections and overlaps. The workow to construct these connections is summarized as follows: (1) Construct initial connections and arrange to avoid overlap.
(2) Optimize the curve order, n, based on an input parameter, s, used to control the channel sparsity.
(3) Optimize the channel width, w, and channel length, h, based on the desired total channel length, L 0 .
(4) Build an AABB for each connection and check for intersections using AABB collision detection algorithms.
Renement of square curves and other interactive functions have been incorporated into AFAN. Microvessels were rendered as truncated generalized cones (TGCs) 61 dened by two adjacent points along a ber that outline the local vessel shape and diameter (Fig. 7d). Volumetric ow and pressure were interactively calculated and visualized using arrow glyphs color-    mapped by ow speed. A group of keyboard and mouse commands are registered for interactive purposes, allowing users to create customized augmented channels attached to an RMN. Conveniently, boundary pressure values are specied with only a few clicks.

Two-dimensional Wheatstone bridge network
We rst designed a simple planar network based on the Wheatstone bridge structure which allows for control of the ow direction in the central channel based on varying boundary conditions. We used AFAN to specify two sets of boundary conditions with inverse ow directions in the central channel via different augmented connections (Fig. 5). We imported the resulting volumes into COMSOL (https://www.comsol.com) to simulate a laminar ow using the same boundary setups. The comparison shows consistency of ow directions between two results, demonstrating our ability to control the ow prole within the ROI network via synthetic connections, or augmentation channels.
Aer fabrication, analysis of FITC-dextran images (Fig. 6a) showed a close agreement in vessel diameters between the mask and fabricated network as previously demonstrated. 31 Segment diameters range from 18.2 mm to 19.2 mm in segments   (Fig. 6c). Segment 2 shows a slightly higher value than the simulation; this is likely due to particle lodging in the bifurcation between segments 1, 2, and 3 that causes a slight narrowing at the inlet of segment 2. Nevertheless, this experiment demonstrates that our AFAN soware can be used to design and model physiologically relevant microuidic networks. Here, we are referring to relevance with respect to wall shear stress that averaged 45 dynes per cm 2 in the synthetic network which closely matches measured in vivo values for capillaries. 55 Three-dimensional, augmented, biomimetic network Microvascular networks are difficult to accurately reconstruct due to the high spatial resolution and large volumes required (Fig. 7b). Microvessels within these networks are oen less than 10 mm in diameter (Fig. 7c), 62 however large volumes are necessary for understanding connectivity patterns and identifying the desired ROIs within microvascular beds. 63 Images were therefore acquired using knife edge scanning microscopy (KESM) 46 at a resolution of 0.6 mm laterally and 1.0 mm axially, which is sufficient to resolve the smallest microvessels. 47 This method provides whole organ images at high contrast, making them easier to reconstruct. Vessel center lines and radii were extracted using an automated predictor-corrector algorithm 49 that reconstructs the medial axis of each capillary ber. This created an explicit graph model storing the architecture of the segmented network (Fig. 7d) and provided an accurate description of connectivity for ow simulations. The average velocity of each channel was initially estimated using AFAN (Fig. 8a) while the maximum velocity of each channel was calculated using SimScale (Fig. 8b). The comparison shows that our linear method succeeds in quantifying the velocity and pressure elds across the augmented network ( Fig. 8c and d). It also shows that the deviation between the average ow velocity and maximum velocity can be compensated by a scaling factor, which is equal to 2 for circular channels (Fig. 8c). This is the result of integrating the velocity prole along the channel cross-section.
An RMN with vessel diameters ranging from 10-20 mm (Fig. 2) was also designed and fabricated. A 3D volume rendering of the fabricated network ( Fig. 9) was created using Amira (ThermoFisher Scientic), which demonstrates a great agreement in network structure between the mask and fabricated network.

Conclusion
In this paper, we proposed a new method for designing microvascular networks to enforce controlled ow when fabricated in microuidic devices. To fulll this aim, we developed soware for simulating, characterizing, and visualizing microvascular networks for in vitro applications. We demonstrated the ow predictability of our linear model using rigorous CFD methods. We also demonstrated the viability of our network augmentation technique by fabricating an AFAN-developed microuidic device and tracking uid ow. The predicted network characteristics, including the structure, uid velocity, and wall shear stress, closely match experimental results. This opens the door to creating microuidic models of microvascular networks whose structure and uid ow parameters mimic that of their in vivo counterparts. This will enhance the use of in vitro cell culture models by allowing researchers to more closely recapitulate in vivo ow patterns. Ideally, we can bridge different sub-networks in series to study blood circulation between different regions or across whole organs, such as the brain. 64,65 AFAN has been integrated into an open-source application available online (https://www.stim.ee.uh.edu). However, this current approach is based on the assumption of Newtonian ow, and thus limited to microvessels with high shear rates. To fully represent realistic in vivo ow, we plan to explore future modications to our algorithm that will account for non-Newtonian effects, as well as validate our results using more accurate ow media.

Authors' contributions
Jiaming Guo implemented the AFAN soware, optimized the network augmenting method, and conducted CFD-based simulations. Pavel Govyadinov segmented and reconstructed the microvascular data. Keely Keller and John Slater fabricated the microuidic networks and conducted the particle tracking experiments. Paul Ruchhoe and David Mayerich developed the theoretical model. All authors reviewed and contributed to the nal manuscript.

Conflicts of interest
There are no conicts to declare.