Rapid model-guided design of organ-scale synthetic vasculature for biomanufacturing

Abstract

Our ability to produce human-scale biomanufactured organs is limited by inadequate vascularization and perfusion. For arbitrarily complex geometries, designing and printing vasculature capable of adequate perfusion poses a major hurdle. We introduce a model-driven design platform that demonstrates rapid synthetic vascular model generation alongside multifidelity computational fluid dynamics simulations and three-dimensional bioprinting. Key algorithmic advances accelerate vascular generation 230-fold and enable application to arbitrarily complex shapes. We demonstrate that organ-scale vascular network models can be generated and used to computationally vascularize >200 engineered and anatomic models. Synthetic vascular perfusion improves cell viability in fabricated living-tissue constructs. This platform enables the rapid, scalable vascular model generation and fluid physics analysis for biomanufactured tissues that are necessary for future scale-up and production.

Fig. 1.. Perfusion performance and capabilities of the synthetic vascular toolkit.

(A and B) Representative lattice and synthetic vascular networks at a fixed total blood volume in convex (cube) and nonconvex (biventricle) domains. (C) Net resistance per unit length (∂R_net) for networks in (i) convex and (ii) nonconvex geometries, normalized by the maximum value (∂R_{net, max}) and plotted against the number of vessels (N_vessels). Lat., lattice; Vas., synthetic vasculature. (D) Minimum flow ranges [(Q_min/Q_inlet), 1.0] versus network size in (i) convex and (ii) nonconvex geometries. (E) Normalized lower quartiles of solute concentrations (C_LQ/C_inlet) within vessels for (i) convex and (ii) nonconvex network geometries. (F) Average bulk solute concentrations in tissue (⟨C_tissue⟩ / C_inlet) approximated by a linear depletion model, plotted for (i) convex and (ii) nonconvex tissues over increasing network sizes (solid lines indicate moving averages). (G) Color-mapped tissue concentrations for networks of three different sizes in (i) convex and (ii) nonconvex geometries, with average bulk concentrations and vessel counts indicated. (H) A biventricle tissue domain (Ω^t) featuring challenging thin-walled regions (ventricular septum, right ventricular wall). (I) Open-loop synthetic vasculature tree with 1000 outlets defining the fluid domain (Ω^f) for blood or media flow. The root (seed) location is indicated. (J) Matched CFD simulation generated from the 1000-outlet tree, with an inlet flow waveform and outlet impedances modeling capillary and venous resistance. (K) Connections between synthetic arterial and venous trees form a closed-loop network suitable for embedded 3D printing.

Fig. 2.. Performance of synthetic vascular acceleration techniques.

(A) Conceptual schematic of partial binding for bifurcation optimization. A “constructor” builds functions for bifurcation ratios (F_β), hydraulic length updates (F_L*), and hydraulic resistance updates (F_R*) based on locally bound tree configurations (dependent on bifurcation position x). The resulting cost function (F_cost) is returned for optimization. (B) Theoretical scaling complexities (red) for previously reported bifurcation optimization schemes versus partial-binding methods (blue). Plots reference bifurcation depth in full, balanced trees (shading: ±2 SD, N = 100). (C) Alignment accuracy versus a high-resolution finite-difference solution for partial binding (blue) and a prior approximation (red). Solid lines show median alignment (shading spans minimum-maximum alignment; N = 100). All trees are generated to 8000 outlets within a cube perfusion volume. (D) Schematic of partial implicit volumes. Surface point data are decomposed into local patches, then blended during implicit reconstruction. (E) Partial implicit volumes generate vascular trees (1000 outlets) across (i) a cube, (ii) an annulus, (iii) a biventricle heart, and (iv) a brain gyrus. (F) Scaling complexity for mesh-based versus implicit evaluations (shading = ±2 SD, N = 25). (G) Collision avoidance schematic. Sphere proximity detects potential collisions for precise OBB checks. (H) Sphere proximity alone is cheap but less accurate (error bars ±2 SD, N = 100). (I) Used as a filter, heuristic sphere proximity eliminates most OBB evaluations as tree complexity grows (error bars ±2 SD, N = 100).

Fig. 3.. Automatic multifidelity hemodynamics modeling.

(A) Major steps for automatic 3D watertight model generation from discrete synthetic vessels. (B) Normalized pressure ⟨P⟩, velocity along streamlines ⟨∣v∣⟩, and wall shear-stress magnitude ⟨∣τ_wall∣⟩ obtained from 3D finite-element method simulation with a steady inflow profile. (C) Illustration of reduced-order 1D and 0D models extracted from original 3D model files. Average normalized pressure ⟨P⟩. (D) and flow rates ⟨Q⟩ (E) obtained from 0D steady-flow simulations for vascular trees with 1000 terminals within (i) cube, (ii) annulus, (iii) biventricle, and (iv) brain gyrus tissue domains.

Fig. 4.. Scalable vessel densities in algorithmic generation and bioprinting of vascular models.

(A and B) Synthetic vascular trees scaled across three terminal vessel densities (10⁴, 10⁵, 10⁶) terminal vessels shown for (A) biventricular and (B) annulus geometries (scale bars are 1 cm). (C and D) (i) Tissue (gray) and vascular domains (red and blue) depicted as lateral and axial views alongside corresponding models with 500 interconnected unevacuated branches, printed in Carbopol ink containing red pigment. (ii) Orthogonal views of printed (C) biventricular and (D) annulus geometries (scale bars are 1 cm).

Fig. 5.. Perfusion and viability of bioprinted vascular models.

(A) Two-step FRESH printing and (B) OCT gauging of printed construct. (C) Dye perfusion demonstrating watertight lumen. (D) CFD simulation (Q_inlet = 0.25 ml/min) showing (i) pressure and (ii) wall shear stress. (E) Schematic of PIV validation setup. (F) Normalized flow rates (q/q_in) from PIV compared across three vessels using unpaired t tests with Welch’s correction (N = 3). Averaged (dashed) and individual (“+”) CFD simulation data are shown. Error bars indicate ± 1 SD. n.s., not significant. (G) 3D vascular model schematic (boundary conditions: inlet Γ_inlet, outlet Γ_outlet; domains: tissue Ω^t, fluid Ω^f). (H) (i) Pulsatile flow simulation boundary conditions (mean Q_inlet = 0.25, amplitude = 0.1 ml/min). (ii) Simulated (dashed) versus theoretical waveforms (solid) are shown. (I) Flow and pressure along cross sections of a vessel. The inset shows the selected vessel, with numbers indicating cross section slice indices. (J) Printed, unperfused network (scale bar is 1 cm). (K) FITC-Dextran perfusion at varying flow rates (scale bars are 1 cm). (L) Volumetric model and circularity after ink evacuation (scale bars are 1 mm). (M) μCT reconstruction and deviation plot (scale bars are 1 cm). (N) Nested bioprinting schematic of perfused annulus. (O) Printed annulus matrices: (left) annulus alone and (right) combined vascular and tissue (scale bars are 1 cm). (P) Vascularized annulus viability: (i) sectioning plane (top) and cultured annulus after 7 days (bottom; scale bar is 1 cm); (ii) LIVE/DEAD (CaAM/EthD-1) staining (top), with normalized radius (R) and core region ∈ [0.2 to 0.8] indicated for average radial cell viability (bottom) (N = 3 tissue constructs per condition, scale bars are 1 mm); (iii) comparison of integrated viable cell count [perfused: (417 ± 37)-fold higher] and viability (perfused: 37.8 ± 33.6%; media control: 0.4 ± 0.4%) within the core region normalized to media control using unpaired t tests with Welch’s correction (N = 3; ** indicates 0.05 < p < 0.01); (iv) Kullback-Leibler divergence comparing observed and approximating distributions; (v) spatial distribution of viable cells near channels; and (vi) polar histogram of viability near the perfused channel (N = 3).