Closed-loop vasculature network design for bioprinting large, solid tissue scaffolds

The vasculature network generation model proposed in this study comprises four major steps—medical image processing, vascular network generation, mass transport modeling and new tether points generation, as shown in figure 1. The process commences with medical image processing and segmentation to obtain the 3D structure of the target tissue or organ. The medical images in the form of a stack of slice images representing the cross-section structure of the tissue/organ at different depths were used in this study. The target tissue/organ regions were segmented in each slice image and stitched to form a point cloud representation of the tissue/organ volume. This point cloud was then translated into a 3D model of the tissue.

Zoom In Zoom Out Reset image size

This study used an open-source dataset of abdomen MRI with the liver as the target organ for designing the vasculature. The MRI dataset was retrieved from The Cancer Imaging Archive (TCIA, www.cancerimagingarchive.net/). Characteristics of the MRI dataset used to generate the liver 3D model are listed in table 1. We agree that variability among datasets can exist due to variations in the MRI scanners, scan parameters, patients, scanned organs and the type of imaging sequence. We believe the selected dataset provides a suitable representative case to demonstrate the proposed design framework. The selected open-source database TCIA is a well-curated database providing essential and reliable resources for medical imaging research. We aimed to demonstrate our method on a thick and dense organ with indispensable vasculature requirement and selected the abdomen MRI dataset owing to its completeness and high quality. Moreover, MR image dataset provides information about tissue water content which can be used to infer various tissue properties. Different MR modes can simplify the segmentation of the target tissue/organ structure from the images. However, the vasculature design method developed in this study is kept generic to allow the use of other imaging modalities. Further, MR elastography data can be used to include tissue viscoelastic properties in guiding the vascular network generation. Different image modes of the MRI dataset images are shown in figure 2(A). T1- and T2-weighted images with and without gadolinium contrast agents were recorded.

Zoom In Zoom Out Reset image size

The MR slice images were segmented based on intensity thresholding to isolate the liver region represented by higher pixel intensity (gray level value). The threshold value was determined manually, however, other improved medical image segmentation methods can be applied. The segmented sections of the MR images were stitched to form a 3D point cloud dataset which was further filtered by applying a point clustering operation. Euclidean distance between the points was used as the clustering criteria. Small sub-clusters and points lying away from the clusters were removed to form an initial 3D liver model. It should be highlighted that the MRI images contain several regions of similar signal intensity that hamper intensity-based segmentation approaches and outputs point cloud data with undesired noise levels. Therefore, manual segmentation of target tissue/organ from the slice images by marking the tissue outlines can be used as an alternative approach. The 3D model achieved from the filtered point cloud was then converted to a 3D solid mesh comprising tetragonal elements. Since the mesh includes several defects that prevent its use in numerical simulations, mesh smoothing and repair steps were applied to improve the mesh quality. The final liver 3D model was then sliced to obtain a stack of binary images showing the cross-section at different depths. The slice spacing parameters specified in the MR image dataset (table 1) were utilized in this step. In this study, the segmented images were manually examined by an expert user, however, to further improve the accuracy and clinical relevance, the segmented images and the obtained organ 3D model can be rechecked by medical experts. The obtained slice images were used as digital masks and superimposed with the corresponding input MR images to extract the MR image data in the superimposed regions. A stack of segmented MR images was obtained with this step to provide gray-level intensity at each voxel of the generated liver 3D mesh. Next, the intensity data from T2 mode images from the segmented slice images was used for estimating a porosity equivalent parameter ${\eta _{{\text{porosity}}}}$. The T2 mode images specifically selected as the pixel gray level intensity in T2 mode are proportional to the tissue water content. Since the designed scaffolds are intended for biofabrication with porous hydrogels, the approximate porosity of gelatin methacrylate (GelMA) hydrogels was used in this study. GelMA hydrogels have been reported to exhibit a porosity in the range of 80% to 90% [36, 37] in the literature which was used to define a linear ${\eta _{{\text{porosity}}}}$ variation as shown in equation (1)

$$\begin{equation}{\eta _{{\text{porosity}}}} = { }0.8 + \left[ {\left( {0.9 - 0.8} \right) \times \left( {\frac{{{I_{x,y,z}}}}{{I_{x,y,z{ }}^{{\text{max}}}}}} \right)} \right],\end{equation} \tag{ 1 }$$

where ${I_{x,y,z}}$ represents the voxel intensity at (x,y,z) coordinates from the liver volume formed by T2 mode images and $I_{x,y,z{ }}^{{\text{max}}}$ represents the volume maximum intensity.

The histogram of the voxel gray level intensity of the liver volume formed by T2 mode images was obtained, corresponding to the distribution of the voxel water content in the liver. A threshold was applied to identify high water content regions (i.e. regions with high ${\eta _{{\text{porosity}}}})$. These regions were then filtered from the liver volume and converted to a point cloud on which a point clustering step was performed to form multiple sub-clusters using Euclidean distance criteria. Centroids of these sub-clusters were identified as an initial set of tether (attraction) points for guiding the vascular network formation. A pair of initiation node points was also identified on the liver 3D model surface to initiate the vascular network formation. The described sequence of steps in the MR image processing with dependencies and outputs is graphically represented in figure 2(B).

The tether points and initiation points extracted from MR image processing were transferred to a custom venation program created in MATLAB^™. Figure 3 shows a graphical representation of the different mathematical steps involved in the venation algorithm. The algorithm commences by computing the distance between the first node (initiation points ${I_1}\,$or ${I_2}$ and the tether points (figure 3(A)). Several tissues and organs can include a cavity or other heterogeneous components (e.g. bones where the designed vasculature cannot be directly printed); therefore, a domain crossing restriction was introduced to restrict the vasculature within the feasible domain of the tissue/organ. To implement this, the lines connecting a node with all available and feasible tether points were formed and examined for intersection with the boundaries of the tissue/organ 3D structure. Any intersection led to the removal of the corresponding tether from the set of feasible tether points for the current node and was considered to not attract the growing vasculature network towards it. This criterion was re-evaluated as the vascular network evolved and presented as a possibility of the network moving around the cavity or bone volumes. A list was populated containing pairs of all the tether points attracting the current node to advance the vasculature. A bi-directional approach was used in determining such pairs where a node was attracted by all feasible tether points, but a tether point was only able to influence the geometrically closest node. The intersection of these sets provided the required tether-node pairs. Using the respective coordinates of the node and the tether points, the unit vectors at the current node were obtained pointing towards the attracting nodes. All the vectors were combined to compute a resultant vector that provided the stepping direction for the growing vasculature and the location of subsequent node. It should be highlighted that the stepping step can be evolved to scale the respective unit vectors and prioritize the attracting tether points based on the tissue physical characteristics as their location. This capability can aid the vasculature design process in various complex tissues where several factors (e.g. presence of different cell types, tissue composition with high fat content) affect the vasculature. The computed resultant stepping vector was then used to place a new node ${n_1}\,$or ${n_2}$. The new node was added to the growing vascular network and the line joining the new node and the previous node (${I_1}\,$or ${I_2}$) integrated with the growing branch as shown in figure 3(B). With the updated vascular network, the pair of feasible attracting tether points for each node were determined again with the methodology described earlier.

Zoom In Zoom Out Reset image size

As the branches grew towards the tethers, the influence of the tether points was limited to the geometrically closest nodes combined with an additional criterion described later. For example, the ${n_1}$ node was closer to several tether point as compared to ${I_1}$, which resulted in omission of those tether points when determining the feasible attracting tether points for node ${I_1}$. This exclusion was determined by examining the relative neighborhood graph of the tether points as explained later in this section. Considering subsequent steps when the network had further evolved and included the nodes ${n_1}$, ${n_2}$, ${n_3}$ and ${n_4}$ as shown in figure 3(C). At this stage, the tether-node pairs were determined for all the nodes and only the tether ${s_p}$ was found to be an eligible attracting tether point for the node ${n_1}$. This resulted in the sprouting of a new branch at the node ${n_1}$ approaching the tether ${s_p}$ and created a bifurcation in the growing vascular network (figure 3(C)). As the vasculature branches get closer to a tether point, a termination criterion was applied such that if the separation between a node and the tether point is less than or equal to the defined kill distance, the tether and the node were joined and terminated. Termination of the tether inactivates its influence over any nodes in further computation and the terminated nodes mark the end of the growing vasculature branch. When the branches from nodes ${n_1}$ and ${n_3}$ approached the tethers ${s_p}$ and ${s_q}$ and entered within the kill distance, the tethers ${s_p}$ and ${s_q}$ and the nodes ${n_1}$ and ${n_3}$ were terminated (figure 3(D)). On the other hand, all the branches approaching a tether point were tracked when a closed vascular network was sought. When a node was placed within the kill distance of the tether point, it was inactivated. However, the tether point was inactivated only when all the approaching branches had been terminated or the approaching branch nodes fell out of the tether point's relative neighborhood and were excluded from its influence in attracting the branch. The geometric interpretation of the various computation steps involved in determining the relative neighborhood of the tether S are shown graphically in figure 3(E). In this illustration, a scenario is considered when four branches approach the tether S along different directions. The leading nodes of these approaching branches were designated as nodes ${n_A}$, ${n_B}$, ${n_C}\,$ and ${n_D}$, respectively. The blue circle marks the locus of a point located at a distance of ${d_{SD}}$ (i.e. the distance between S and ${n_D}$) is indicated by a blue circle. When the leading nodes approaching the tether lie within the same quadrant, the node geometrically closer to the tether lies in the relative neighborhood of the tether and the other nodes are omitted. However, when the leading nodes on the approaching branches lie on different quadrants, the geometrical proximity is not the sole criteria for determining the nodes in the relative neighborhood of the tether. In such a scenario, a node ${n_i}$ was considered to be in the relative neighborhood of the tether S if the separation between the tether S and any other node ${n_j}$ was less than the separation between S and ${n_i}$ but the separation between ${n_i}$ and S was less than the separation between ${n_i}$ and ${n_j}$ [34]. A generalized form of this criterion is mathematically represented by equation (2).

$$\begin{equation}\left( {\forall {n_j} \epsilon N} \right)\|{n_i} - S\| < {\text{max}}\left\{ \|{{n_j} - S\|,\|{n_i} - {n_j}}\| \right\}\end{equation} \tag{ 2 }$$

After applying this definition of relative neighborhood, the nodes ${n_A}$, ${n_B}$ and ${n_D}$ were included in the relative neighborhood of the tether S while the node ${n_C}$ was omitted because ${d_{SC}} > {\text{max}}\left\{ {{d_{SD}},{d_{DC}}} \right\}$. Intrinsically, this criterion forced a geometric constraint that multiple branches in the vasculature can simultaneously approach a tether only if the lines joining the leading nodes of these branches and the approached tether point are at obtuse angles with respect to each other (i.e. the angle between two branches approaching the same tether cannot be less than $\left( {\frac{\pi }{2}} \right){\text{rad}}$). This geometric constraint also imposes a limit on the vasculature density. Therefore, we revised this constraint in order to generate denser vascular networks and equation (2) was improved to include a parameter ${\varphi _{{\text{approach}}}}$ (equation (3)). For instance, if $\theta$ is set to $\left( {\frac{\pi }{3}} \right){\text{rad}}$, multiple vasculature branches could simultaneously approach the tether point with a minimum angle of $\left( {\frac{\pi }{3}} \right){\text{rad}}$ between them and output a denser vascular network. On the other hand, increasing $\theta$ beyond $\left( {\frac{\pi }{2}} \right){\text{rad}}$ will result in the formation of a sparser vasculature.

$$\begin{align}&\left( {\forall {n_j} \epsilon N} \right)\|{n_i} - S\|\nonumber\\&\quad < \max \left\{ \|{{n_j} - S\|,\left( {{\varphi _{{\text{approach}}}} \times \|{n_i} - {n_j}}\| \right)} \right\}{ }\end{align} \tag{ 3 }$$

where ${ }{\varphi _{{\text{approach}}}} = \left( {\frac{\pi }{{2\theta }}} \right)$

Once the vascular network was generated from the venation algorithm, a cleaning operation was performed to remove short branches that included less than two nodes and did not connect to a tether. While designing a closed vascular network, the cleaning operation identified and removed all the open branches (i.e. not connected to other branches directly or at tethers) in the vascular network. The tethers that were connected to a single branch were also removed in this process along with the branch. As the vasculature comprised of several splits where the parent branch divided into multiple child branches, the Murray's Law with $n = 3$ was applied to determine the thickness of the respective branches (equation (4)).

$$\begin{equation}r_{{\text{parent}}}^n = {\mathop \sum \nolimits }r_{{\text{chil}}{{\text{d}}_{\text{i}}}}^n\end{equation} \tag{ 4 }$$

The thickness of the primary branch of the vasculature at the initiation nodes was defined manually, and the thickness of subsequent child branches was computed with equation (4). When determining the branch thickness in closed vascular networks, Murray's Law poses an ill-defined problem with no solution. Therefore, an additional constraint was applied by restricting the minimum branch thickness to the dimensions of the smallest features that can be fabricated by the bioprinting system that will be used subsequently for fabricating the designed vascular tissues. The entire vasculature was a connected network of nodes which were points in a 3D cartesian space. The branch thickness was applied by placing concentric circles of determined diameter at and between the connected nodes. After applying the branch thickness to the entire vasculature, a point cloud was formed that was transformed into a 3D mesh object to be used in subsequent steps or for biofabrication.

The obtained vasculature point clouds for both artery and vein were imported into MeshLab, an open-source meshing software. Point normals were computed for the point clouds, followed by screened Poisson surface reconstruction to generate watertight surfaces. Parameters such as reconstruction depth, interpolation weight and Gauss Seidel relaxations were adjusted to create a smooth surface for the artery and vein 3D models. The resulting surface mesh was imported into Autodesk Meshmixer and converted into a solid tetragonal mesh. Artery and vein mesh, when assembled, were found to intersect at a few points that were readjusted to create a combined intersection free mesh. As these meshes consist of approximately 50,000 faces, intersecting faces and inward point facing normals occur frequently resulting in mesh with poor quality. Thus, several smoothing and remeshing operations are recommended before proceeding with further simulation steps.

Current studies model and fabricate networks that have arteries and veins connected together (i.e. closed vasculature) without the presence of capillaries, which creates a shunted blood flow [18, 19, 21, 38]. Although these veins and arteries pass through the fabricated tissue, the diffusion rate of nutrients from the walls of such networks is significantly slower than the fluid flow within the network. In the present model, the veins and arteries are not directly connected (i.e. open vasculature). Instead, nutrient transport occurs through the porous hydrogel between the arteries and veins. Nutrient transport through these pores crudely mimics a cluster of tiny capillaries. The described venation algorithm was tested to design a closed shunted vasculature forming around the brain, as shown in figure S1(A). The capability to design an open vascular architecture was also demonstrated by designing a combination of arterial and venous networks inside a liver (figure S1(B)).

The processed 3D mesh combining the artery and vein structure and the 3D model of the liver generated from the MR image processing was imported into COMSOL to simulate the mass transport of various cell growth media components. We focused on the mass transport characteristics of six components—glucose, glycine, glutamine, riboflavin, human serum albumin (HSA) and oxygen. These components represented the classes of sugar, amino acid, vitamin, protein, and gas, respectively. The designed 3D liver structure and the vascular network are intended to be bioprinted with porous hydrogels. Therefore, the surface of the artery and vein network within the liver was defined as an internally permeable wall allowing the diffusion of cell growth media components. The vascular networks were defined as hollow domains filled with cell growth media. The vascular network's primary branches beginning at the initiation nodes are, by design, located within the tissue domain, which makes them inaccessible to apply required boundary conditions. Therefore, the primary initiation branches were extended to protrude out of the liver volume. The extended protruded parts in the arterial and venous vascular networks were defined as hollow tubular structures for the inlet and outlet of the cell growth media (figure 4(A)). Next, we defined the liver body as a porous medium saturated with an aqueous solution, and the initial concentration of various cell growth media components was set as zero. Since a porous hydrogel will be utilized to fabricate the liver scaffold, it is expected to allow diffusion of various chemical species through the outer surface. However, an epithelial cell layer forms the outer surface of the organs and prevents such diffusive mass transport and limits the mass transport contribution to the vasculature. Therefore, we modeled the external wall of the liver body to be surrounded by mucous, which prevents a diffusive mass transport by applying a no flux boundary condition. The local porosity ${\epsilon _P}$ within the liver body was defined using the porosity equivalent parameter ${\eta _{{\text{porosity}}}}$ determined earlier from the MR images. The spatial porosity variation in y-z and x-y cross-sectional planes is shown in figures 4(B)–(D). The diffusive mass transport in the liver domain was modelled using equation (5) by considering a saturated porous media.

$$\begin{equation}\frac{{\partial \left( {{\epsilon _P}{c_i}} \right)}}{{\partial t}} + \nabla \cdot {J_i} = { }{S_i},\,{\text{where }}{\epsilon _P} = {\eta _{{\text{porosity}}}}\end{equation} \tag{ 5 }$$

The diffusive flux ${J_i}$ was defined by equation (6) using an effective diffusivity ${D_{e,i}}$, considering a negligible contribution in mass transport from dispersion.

Zoom In Zoom Out Reset image size

$$\begin{equation}{J_i} = \left[ {\left( {{D_{e,i}}} \right)\nabla {c_i}} \right]\end{equation} \tag{ 6 }$$

The effective diffusivity ${D_{e,i}}$, in turn, was computed using equation (7) by using Millington and Quirk model. This model defines the fluid tortuosity ${\tau _L}$ as a function of media porosity ${\eta _{{\text{porosity}}}}$, and ${D_{L,i}}$ represents the diffusion coefficient of the various cell growth media components in the aqueous medium

$$\begin{equation}{D_{e,i}} = \frac{{{\eta _{{\text{porosity}}}}}}{{{\tau _L}}}{D_{L,i}},{\text{ where}}\,{\tau _L} = \eta _{{\text{porosity}}}^{-\frac{{\text{1}}}{{\text{3}}}}\end{equation} \tag{ 7 }$$

The thin permeable vasculature walls allowed a diffusive mass transport with nutrient diffusion flux ${J_{{\text{wall}}}}$ given by equation (8), where the concentration of various components in the liver domain and the vasculature domain are denoted by ${[{c_i}]_{{\text{liver}}}}$ and ${\left[ {{c_i}} \right]_{{\text{vasculature}}}}$, respectively, and ${\delta _{{\text{wall}}}}$ represents the thickness of the thin permeable outer wall of the vasculature.

$$\begin{equation}{J_{{\text{wall}}}} = {D_{L,i}} \times \left( {\frac{{{{\left[ {{c_i}} \right]}_{{\text{liver}}}} - {{\left[ {{c_i}} \right]}_{{\text{vasculature}}}}}}{{{\delta _{{\text{wall}}}}}}} \right)\end{equation} \tag{ 8 }$$

Convective mass transport was not included in this model as the growth media flows through the vasculature and the diffusion of various chemical species occurs through the permeable wall of the vasculature channels and beyond the vasculature. The species transport is not influenced by the growth media flow. Although the cell growth medium comprises numerous components, for simplicity, five components were considered, and their respective concentrations and diffusion coefficients were obtained from previous experimental investigations (table 2). The entire liver scaffold was considered to encapsulate the hepatocytes as a bioprinted construct. As the cell media components diffused through the liver, they were also consumed by the encapsulated cells with the consumption rates listed in table 2. The local consumption rate ${S_i}$ of these components was computed by equation (9),

$$\begin{equation}{S_i} = { }{S_{{c_i}}} \times {\rho _{{\text{cell}}}}\end{equation} \tag{ 9 }$$

Table 2. Diffusive mass transport modelling parameters.

Parameter	Value	Description	References
${D_{{\text{Glucose}}}}$	$0.96 \times {10^{ - 9}}\left( {{{\text{m}}^{\text{2}}}{\text{s}}^{-1}} \right)$	Diffusion coefficient of glucose in water	[39]
${D_{{\text{Glycine}}}}$	$0.99 \times {10^{ - 9}}{ }\left( {{{\text{m}}^{\text{2}}}{\text{s}}^{-1}} \right)$	Diffusion coefficient of glycine in water	[40, 41]
${D_{{\text{Glutamine}}}}$	$1.0 \times {10^{ - 9}}{ }\left( {{{\text{m}}^{\text{2}}}{\text{s}}^{-1}} \right)$	Diffusion coefficient of glutamine in water	[42]
${D_{{\text{Riboflavin}}}}$	$0.42 \times {10^{ - 9}}{ }\left( {{{\text{m}}^{\text{2}}}{\text{s}}^{-1}} \right)$	Diffusion coefficient of riboflavin in water	[43]
${D_{{\text{HSA}}}}$	$0.06 \times {10^{ - 9}}{ }\left( {{{\text{m}}^{\text{2}}}{\text{s}}^{-1}} \right)$	Diffusion coefficient of human serum albumin in water	[44]
${D_{{O_2}}}$	$2.69 \times {10^{ - 9}}\left( {{{\text{m}}^{\text{2}}}{\text{s}}^{-1}} \right)$	Diffusion coefficient of oxygen in water	[45]
${\left[ {{C_{{\text{Glucose}}}}} \right]_o}$	${\text{24}}{\text{.9 mM}}$	Initial concentration of glucose in growth medium	[46, 47]
${\left[ {{C_{{\text{Glycine}}}}} \right]_o}$	${\text{0}}{\text{.4 mM}}$	Initial concentration of glycine in growth medium	[46, 47]
${\left[ {{C_{{\text{Glutamine}}}}} \right]_o}$	$0.4{\text{ mM}}$	Initial concentration of glutamine in growth medium	[46, 47]
${\left[ {{C_{{\text{Riboflavin}}}}} \right]_o}$	$0.001{\text{ mM}}$	Initial concentration of riboflavin in growth medium	[46, 47]
${\left[ {{C_{{\text{HSA}}}}} \right]_o}$	$0.0038{\text{ mM}}$	Initial concentration of human serum albumin in growth medium	[48]
${\left[ {{C_{{O_2}}}} \right]_o}$	${\text{0}}{\text{.2 mM}}$	Initial concentration of oxygen in growth medium	[49, 50]
${S_{{\text{Glucose}}}}$	$3.34 \times {10^{ - 18}}\left( {{\text{mol/cell}}}\times {\text{s}}\right)$	Estimated glucose consumption rate	[51]
${S_{{\text{Glycine}}}}$	$3.93 \times {10^{ - 18}}\left( {{\text{mol/cell}}}\times{\text{s}}\right)$	Estimated glycine consumption rate	[52]
${S_{{\text{Glutamine}}}}$	$0.93 \times {10^{ - 17}}\left( {{\text{mol/cell}}}\times {\text{s}} \right)$	Estimated glutamine consumption rate	[53]
${S_{{\text{Riboflavin}}}}$	$5.5 \times {10^{ - 21}}\left( {{\text{mol/cell}}}\times {\text{s}} \right)$	Estimated riboflavin consumption rate of glucose	[54]
${S_{{\text{HSA}}}}$	$3.63 \times {10^{ - 20}}\left( {{\text{mol/cell}}}\times{\text{s}} \right)$	Estimated human serum albumin consumption rate	[55]
${S_{{O_2}}}$	$22 \times {10^{ - 18}}{ }\left( {{\text{mol/cell}}}\times{\text{s}}\right)$	Estimated oxygen consumption rate	[56]
${\rho _{{\text{cell}}}}$	$8 \times {10^6}{ }\left( {{\text{cells/mL}}} \right)$	Initial cell encapsulation density in scaffold	[15, 57]

where ${S_{{c_i}}}$ denotes the molar consumption rate of ${c_i}$ by a single cell per second and ${\rho _{{\text{cell}}}}$ represents the population density of the cells in the scaffold. The cell density was kept as $8 \times {10^6}\left( {\frac{{{\text{cells}}}}{{{\text{mL}}}}} \right)$ based on previous bioprinting studies [15]. The glycine and albumin secretion as a result of the cell metabolic processes was not included in this study for simplification.

The diffusive mass transport in the liver scaffold was simulated up to a duration of 300 min (5 h). Although the cells can survive starvation periods for up to several hours, an upper bound for the starvation period is required in the vasculature design process to constrain the complexity of the vascular networks. This upper bound should be defined by considering the target tissue/organ and the metabolic activity of the constituent cells. In this study, the maximum starvation period was defined as 60 min to impose a stringent condition to drive the branching of the vascular network. The regions in the liver domain without a positive concentration of the cell growth media components at the end of the starvation period were designated as 'starvation' regions. A point cloud dataset was populated with a list of coordinates representing the 3D liver volume and the corresponding concentration of cell growth media components. The starvation regions were then filtered from this dataset, and a point clustering operation with Euclidean distance between the points as the clustering criteria was applied to segment the filtered point cloud dataset. The centroids of the segment clustered were selected as the new tethers for venation. The final list of tethers included the new tethers as well as the previous set of tethers used to generate the vasculature. New tethers in close proximity to the old tethers were removed, and the set of updated tethers was provided as input for the next iteration of the venation to improve the vascular network.

The liver 3D model was obtained from an open-source MR image dataset by applying the methodology described above (section 2.1). The MR images were first segmented by applying an intensity threshold on the T1-contrast images (figure 5(A)) to remove dark regions surrounding the liver tissue. The segmented slice images were then individually converted to a 2D point cloud dataset and were further refined by applying a point clustering operation. Post filtering, the slice images were stacked to obtain a 3D point cloud by applying the voxel scaling parameters from table 1 and then the point cloud was converted to a 3D mesh by applying triangulation operations, as shown in figure 5(D). The generated 3D liver model had a rugged topography with sharp vertices and disconnected faces due to the low resolution across the MRI dataset slices and error in the segmentation of exact liver boundaries. Mesh smoothing and mesh repair operations were performed on the 3D liver model to obtain a smooth surface topography, as shown in figure 5(E). The generated model included the aorta along with the liver tissue, which was removed manually to obtain the final liver 3D model (figure 5(F)). Mesh smoothing and repair operations resulted in slight deformation in the overall 3D liver model and introduced a volume shrinkage of 6% when compared to the rugged topography 3D model shown in figure 5(D). The smoothed liver 3D model was sliced again with the same slicing parameters as the input MRI dataset to generate the binary cross-section images at different heights (figure 5(B)). The corresponding slice images of the liver 3D models before and after the mesh smoothing and repair operations (figures 5(D) and (E)) were compared to compute the reconstruction error introduced in each slice due to various mesh manipulations (figure 5(G)). The cross-section area reconstruction error was maintained below 10% in the majority of the slices. However, the liver cross-section in the inferior slices (i.e. slices at the bottom part of the liver) suffered comparatively large distortion (15%–18%). As this study emphasizes the development of a process workflow, a simple segmentation and 3D organ reconstruction approach is used. However, the accuracy of the image segmentation and organ reconstruction can be improved using previously reported methods along with the inclusion of machine learning models [8–10, 58].

Zoom In Zoom Out Reset image size

The slice images were superimposed from the smooth 3D liver model without aorta on all the image modes in the input MRI dataset (i.e. T1, T2, T1-enhanced and T2-enhanced) and segmented the superimposed regions, as shown in figure 5(C). These images included the intensity variation at each voxel of the liver volume which was subsequently used to evaluate the local tissue's physical properties. figures 6(A)–(D) show 3D volume plots demonstrating the image signal intensity variation in the liver volume for different MRI modes (T1, T2, T1-enhanced and T2-enhanced). The maximum and minimum grayscale intensity values in the 8-bit input MR images correspond to 255 and 0 and are represented by red and blue, respectively, in a rainbow color map. Further, the histogram for these volumetric data was generated to show the distribution of the faction of voxels at different grayscale intensities, as shown in figure 6(E). The T1-contrast and T2-contrast mode images with gadolinium contrast agent provided a nearly overlapping voxel intensity distribution; however, they had subtle differences. This can also be observed when comparing the T1-contrast and T2-contrast enhanced slice images (figure 2(B)), where the intensity of the liver tissue is comparable. However, a clear difference in the intensity of the vasculature can be observed. The T1 and T2 mode images, on the other hand, show a clearly distinguishable voxel intensity distribution. The narrow voxel intensity distribution in the 3D volume formed by the respective image modes implied relatively uniform tissue characteristics in the liver with slight variations. It should be noted that some tissues/organs (e.g. Kidney) have highly heterogeneous architecture and tissue characteristics and need to be approached with a refined strategy for extracting tissue/organ structure and tissue properties. The 3D point cloud formed by the segmented T2 mode images was further used to compute the porosity equivalent parameter ${\eta _{{\text{porosity}}}}$ (figures 4(B)–(D)) at each voxel, as well as to identify the initial tether points for initiating the venation algorithm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A combination of two distinct vascular networks (arterial and venous) was used to represent the overall liver vasculature, and the porous hydrogel matrix between the two networks mimicked a capillary bed as described above (section 2.2). Using the 3D point cloud volume data of ${\eta _{{\text{porosity}}}}$ derived from the T2 images, two distinct sets of tether points were obtained for generating the arterial and venous vascular networks, respectively, as described previously. The initiation tether points for arterial and venous networks were also defined separately and coincided with the position of the primary artery and vein near the surface. Figure 7(A) shows the obtained initial set of tether points for the venous network generation along with the lower half of the liver body for reference. Figure 7(B) shows a partially formed vein network following the venation computation steps described earlier. Once the arterial and venous vascular networks were generated, they were processed to remove very short branches (i.e. comprising less than two nodes) and apply the respective branch thickness at each node computed using the Murray's Law. Upon applying the branch thickness, the finalized vascular networks were transformed to 3D point clouds and then converted to 3D mesh objects. The generated 3D mesh objects of arterial and venous networks are collectively shown in figure 7(C). The arterial network is shown in red color and the venous network is shown in blue color for clear demarcation. The 3D meshes were then combined with the liver 3D model and a rendered image mimicking the physical appearance of a 3D bioprinted liver was generated utilizing the optical properties of a typical hydrogel (e.g. GelMA) as shown in figure 7(D).

The 3D models of the liver, artery and vein network were imported and combined in COMSOL, and the diffusive mass transport was simulated using the methodology described earlier. Figure 8 shows the isosurface plots representing the concentration of glucose diffused within the liver body from the vasculatures designed in Iteration 1 and Iteration 2. To demonstrate the availability of glucose for the cells and the starvation regions within the liver body, representative isosurface plots at 6, 30, 60 and 300 min, respectively, are shown. A concentration gradient of the diffused species also exists within the isosurface envelope (figure 8 inset image). The extent of diffusion for other species in the cell growth media—glycine, glutamine, riboflavin, human serum albumin and oxygen are also shown in figures S2–S6. The portion of the liver volume covered with the diffused cell growth media components is shown in figure 9(A). The glucose diffusing from the vasculature designed in Iteration 1 covered 56.7% of the liver volume at 5 h while including cellular consumption. On the other hand, 86.8% of liver volume was covered by glucose diffusion in absence of cellular consumption.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Conversely, the vasculature designed in Iteration 2 exhibited improved coverage due to a denser network and covered 64.4% of the liver volume at the end of 5 h in presence of cellular consumption (table S1). However, in absence of cellular consumption, the glucose diffusion from the vasculature designed in Iteration 2 covered 86.5% of liver volume. Further iterations are expected to increase the complexity of the vasculature and achieve complete coverage of the liver volume in the designated maximum starvation duration of 5 h. As various chemical species become available to the encapsulated cells, they begin the uptake of the respective chemical species; however, this uptake is not consistent throughout the liver scaffold and is adaptive in nature, relying on the distribution of the species. The presented model accounts for such adaptive cellular consumption of the growth media components, as demonstrated by figure 9(B), which shows the average consumption rate of various species within the liver scaffold. HSA and riboflavin are present in the cell culture media at very low concentrations and have very low diffusivity in an aqueous medium due to their large molecular size. Concurrently, their cellular consumption rates are very low as compared to other chemical species that result in a similar portion of the liver being covered by albumin and riboflavin as compared to smaller molecules like glutamine and glycine. The adaptive nature of species consumption by the cells, in turn, contributes to the diffusive transport of various species that cause fluctuations in the total portion of liver volume covered (figure 9(A)) in the initial period. We also examined the diffusive transport of various chemical species in the absence of cellular consumption. As shown in figure S7, the vasculature designed in Iteration 2 was suitable to cover a major portion of the liver scaffold with glucose for 5 h.

An often-overlooked pitfall is the differential diffusion of various cell growth media components within the bioprinted constructs. Such differences are neglected in small scaffolds and current studies in which the bioprinted scaffolds are submerged in cell culture media. However, the differences in the diffusivity of various growth media components can lead to significant differences in the composition of unit volume at different regions in large scaffolds. This will be further convoluted by implementing a more physiologically relevant constraint of relying solely on the vasculature rather than submerging the entire scaffold in the cell growth media. This is clearly highlighted in figure 9(A), where the differences in the concentration and molecular size of various species in the cell growth media resulted in striking differences in the coverage portion of the liver body, thus, altering the chemical composition of each voxel.

Although the hepatocytes can sustain starvation periods extending beyond 5 h, we imposed a stricter criterion for increasing the complexity of the designed vasculature in Iteration 2 and beyond by considering the glucose availability at 60 min (figure 10(A)). Such a criterion identifies a larger portion of the liver body as the starved region. The concentration distribution of glucose within the liver body was used to identify the starvation region using a custom MATLAB script and is represented as a point cloud in figure 10(B). These point clouds were further processed by applying a Euclidean distance-based point clustering operation as described in earlier sections. The centroids of these point clusters were then identified and compared with the previous set of vasculature tether points to determine the new set of tether points for designing a vasculature evolved from the previous iteration. The set of tether points obtained for arterial (red) and venous (green) network generation at the end of Iterations 1 and 2 are shown in figures 10(Ci) and (Cii), respectively. The second iteration of the venation algorithm generated a denser network, which in turn led to broader nutrient coverage in the liver body, as shown in figure 8.

Zoom In Zoom Out Reset image size

A vasculature network architecture ideally consists of branched arterial and venous networks and these two networks are connected by a dense capillary bed. Such structure is capable of perfusion under a pressure gradient to generate bulk or interstitial flow of cell growth media. Although arterial and venous networks can be designed, the design of capillary bed poses a challenge. This challenge arises from the different length scales of arteries/veins and capillaries as well as from determining the spatial location of the capillaries. Inclusion of such capillary designs can be achieved by performing multiscale modeling at significantly higher computational complexity. On the other hand, the absence of the capillaries limits the capability to model convective mass transport. In this study, a simplification has been introduced considering that if a bioprinted scaffold is perfused with the cell growth medium, the fluid flows through the connected vasculature and the vasculature walls restrict the flow beyond the vasculature. This leads to a diffusive mass transport occurring from the vasculature walls to the bulk tissue whereas the convective mass transport within the vasculature is significantly larger. Therefore, the cell growth media at any unit volume of the vasculature is replenished quickly and appears to be a constant composition at the diffusive mass transport timescale. This has been modelled in the present study by considering a constant composition of growth media in the vasculature branches and simulating the diffusive mass transport within the bulk tissue. Further, the bulk tissue is modelled as a porous material in which the interconnected pores approximate the presence of a capillary bed. However, including a multiscale model of capillaries can further improve the presented model and provide better insights into angiogenesis and tissue maturation. Moreover, the designed vasculature with simplified architecture is intended to ensure adequate nutrient availability to the encapsulated cells in the entire scaffold in the early stages of maturation after bioprinting. As the vascular scaffold matures, the formation of capillaries will lead to a perfusable vascular network.

The vascular networks generated from the three iterations in this study are shown in figure 11. The perspective views of the network 3D models are shown in figures 11(Ai), (Bi) and (Ci), top views are shown in, (ii), (Bii) and (Cii) and the rendered 3D models are shown in, (iii), (Biii) and (Ciii), respectively. Typically, 4-5 iterations are expected to achieve a volume coverage greater than 80%. Although the presented model captures the consumption of various chemical species, modelling the production of various metabolites by the cells can improve this model. Adaptive cellular consumption rates have been modelled in the presented system; however, such cellular consumption processes can be coupled with known metabolic pathways to define how the cellular consumption kinetics will vary. Moreover, the growth of cells and migration using agent-based models [59, 60] and angiogenesis [61] can also be implemented for a more realistic perspective on the maturation of the fabricated vascularized tissues. The migration and growth of the cells within the scaffold will be governed by the degradation of the biomaterial scaffold which can be implemented in the presented model to improve it. The local physical properties and their variation in the hydrogel scaffold can be predicted based on the bioprinting conditions [62] which, in turn, can be used to predict scaffold degradation as the scaffold matures under the physiological conditions [63] in the extension of the presented study.

Zoom In Zoom Out Reset image size

The scaffold remodeling by encapsulated cells during maturation should include various cell-biomaterial interactions as well as multiple cell types [64, 65]. A suggested pathway can include modeling the enzyme and extracellular matrix secretion by the cells followed by enzymatic degradation of the scaffold hydrogel [66, 67]. These processes can lead to changes in the hydrogel porosity which can be integrated with an angiogenesis model in simulation the diffusive and convective mass transport from a capillary bed between the arteries and veins. Concurrently, the designed vasculature should be experimentally validated. The fabrication of scaffolds with high structural complexity has been demonstrated by evolved biofabrication methods such as embedded 3D bioprinting of hydrogels in a biocompatible gel bath with multiple bioinks, high-resolution stereolithography, or volumetric bioprinting [12, 20, 68, 69]. However, these biofabrication approaches still have specific limitations which can limit their suitability in fabricating the vascular structures designed in this study. We propose a hybrid bioprinting approach combining embedded 3D bioprinting of hydrogels and volumetric bioprinting can potentially achieve the biofabrication of the designed vascular scaffolds. In such hybrid bioprinting approach, the embedded 3D bioprinting of hydrogels can be used to fabricate the vascular network followed by the surrounding tissue fabrication by volumetric bioprinting. For detailed in vitro characterization, small vascular tissues with comparable anatomical complexity can be fabricated with the hybrid bioprinting method. These scaffolds can be subjected to similar physiological conditions in the numerical modeling and the mass transport of various growth media components and local cell growth can be experimentally measured and compared with the simulated values.