A Computational Template for Three-Dimensional Modeling of the Vascular Scaffold of the Human Thyroid Gland

Abstract

We recently designed an innovative scaffold-bioreactor unit for the bioengineering of a three-dimensional (3D) bioartificial human thyroid gland or its miniaturized replica as a part of a microfluidic chip test system. This device is based on the evidence that the 3D geometry of the intraglandular stromal/vascular scaffold (SVS; i.e., the fibrous and vascular matrix) of mammalian viscera plays a key role in guiding growth and differentiation of in vitro seeded cells. Therefore, we initiated a research program focused on computer-aided reconstruction of the 2nd to 4th order intralobar arterial network (IAN) of the human thyroid gland as a reliable surrogate for its 3D SVS, to be used as an input for rapid prototyping of a biomaterial replica. To this end, we developed a computational template that works within the Mathematica environment, giving rise to a quasi-fractal growth of the IAN distribution, constrained within an approximation of the thyroid lobe shape as a closed surface. Starting from edge detection of planar images of real human thyroid lobes acquired by in vivo real-time ultrasonography, we performed data approximation of the lobar profiles based on splines and Bezier curves, providing 3D lobar shapes as geometric boundaries for vessel growth by a diffusion-limited aggregation model. Our numerical procedures allowed for a robust connection between development of lobar arterial trees and thyroid lobe shape, led to a vascular self-similarity consistent with that of a cadaveric lobar arterial cast, and reproduced arterial vessels in a proportion not statistically different from that described for the real human thyroid gland. We conclude that our algorithmic template offers a reliable reproduction of the extremely complex IAN of the adult human thyroid lobe, potentially useful as a computational guidance for bioprinting of thyroid lobe matrix replicas. In addition, due to the simplicity and limited number of morphometrical parameters required by our system, we predict its application to the design of a number of patient-tailored human bioartificial organs and organs-on-chip, including parenchymal viscera and bones.

Impact statement

The study introduces the computer simulation of the three-dimensional (3D) intrinsic vascular matrix of the human thyroid gland, offering a general concept applicable to a number of other human viscera. Indeed, it provides a flexible software tool for reproduction of a 3D surrogate of the organ's 3D stromal matrix, suitable for eventual 3D bioprinting with biomaterials, and recellularization with organ-specific stem cells/progenitors. The final expectation is the design of patient-tailored 3D organ's matrices upon clinical request.

Introduction

It is increasingly recognized that the natural organ's matrix or stromal/vascular scaffold (SVS) plays a critical regulatory role in developmental morphogenesis and related functional maturation of mammalian viscera.¹ In epithelial organs, three-dimensional (3D) tissue patterning, cellular differentiation, vascular morphology, organ size, and hormonal secretory activity are topobiologically influenced by the interactions occurring between SVS and growing cells,² including the thyroid gland.^3,4

Accordingly, we proposed that the 3D geometry of the natural SVS may act per se as a pivotal mechanochemical input for growth and differentiation of thyroid stem cells/progenitor cells.^5–9 Consistently, we showed that a natural SVS obtained by decellularization of the rat thyroid lobe drives in vitro self-assembly and functional differentiation of colonizing rat thyrocytes.^6,9–13

Based on animal models,^7,14,15 recent investigations tackled the in vitro 3D reconstruction of human thyroid tissue. Investigations relied on the capacity of normal and cancerous human thyroid follicular cells and thyroid-differentiated human embryonic and induced pluripotent stem cells to self-assemble into 3D organoids, either spontaneously or driven by a biocompatible 3D scaffold, including a decellularized xenogenic thyroid lobe SVS or exploiting microgravity.^16–20 However, none of these approaches allowed for a robust 3D anatomical reconstruction of a patient-tailored human thyroid lobe.

Therefore, we designed a cell culture bioreactor for self-assembly of a human bioartificial thyroid gland directly on the laboratory bench (i.e., ex situ), including its miniature replica as a part of a microfluidic chip test system. This device includes a core that replicates the 3D architecture of the adult human thyroid SVS fabricated with biomaterials. A basic assumption is that the 3D organ-like (i.e., organomorphic) geometry of the SVS may provide vascular conduits able to supply trophic media and O₂ inside the mass of stem cells/progenitor cells eventually seeded, thus favoring their growth and differentiation up to formation of a complete viable human thyroid gland organoid ex situ.^5,6,9

For this purpose, we initiated a computational and reverse engineering program to identify the best approaches for both replication of the 3D geometry of the adult thyroid SVS and its physical reconstruction by 3D bioprinting/additive layer manufacturing and biomaterials.^21–27 In this study, we present an original core of numerical and rendering routines aimed at virtually reconstructing the intralobar arterial network (IAN) of the human thyroid, chosen as a reliable replica of its entire SVS, constrained within the limits of a thyroid lobe simulation. Preliminary analyses of few building blocks of this computational structure were presented elsewhere.^21,22,25,27

Materials and Methods

Human subjects, 2D imaging of the thyroid lobe, and its 3D simulation

In vivo ultrasonographic images of the human thyroid lobe were obtained from three euthyroid control subjects (a male, two females, age range 17–43 years, body mass index [BMI] range 22.4–24.1 kg/m²) without visible and/or palpable thyroid nodules, screened in one outpatient endocrine clinic (OSTEONET Unit, Medical Center Galliera/CMG—DIMEC, UNIPR, Italy) for a routine examination. Data were collected after approval by the ethical committee of the clinical center, and used under patient's informed consent. A Toshiba Nemium Premio real-time ultrasound machine with 9–12 MHz linear probe was used.

To obtain an approximate representation of the lobe profile, a single lobe projection was chosen, specifically a bidimensional longitudinal section taken at the longest sagittal axis of the lobe. In this manner, we avoided complex computations based on a set of lobe slices, as currently occurring in 3D reconstruction of many human organs by medical rendering techniques. In addition, our approach allowed for determination of features of approximated symmetry (here coined as pseudosymmetry) of the lobe profile, with the aim of obtaining an adequate sampling of points on the lobe contour.

Finally, an analytical description of the entire simulated lobe surface was reached within the computer algebra environment of Mathematica,²⁸ and based on the following computational steps:

Edge detection and spline interpolation of the thyroid lobe profile,

homotopy-like parametrization and 3D reconstruction of the thyroid lobe, and

estimate of the thyroid lobe volume.

Details on each of these steps are available in Supplementary Material and Supplementary Figure S1.

Fractal growth of the intralobar vascular tree through diffusion-limited aggregation

To achieve a 3D reconstruction of the lobar SVS, we focused on the simulation of the intraglandular arterial branches, chosen as a reliable geometric replica of the thyroid SVS. To this aim, we developed an original quasi-fractal algorithm based on the mathematical principle of diffusion-limited aggregation (DLA), leading to a fractal growth as a binary Brownian tree.²⁹ Details are given in Supplementary Material under Auxiliary DLA Subroutine and Template of DLA Subroutines sections including Supplementary Figures S2 and S3.

In addition, to ensure growth of the IAN as a reliable replica of its entire SVS, constraints were imposed on the DLA to ensure that the simulated intralobar vascular tree remained within the boundaries of the reconstructed thyroid lobe, and according to specific tree branching rules. Details on this can be found in Supplementary Material under Constraints on DLA Growth section including Supplementary Figure S4.

Finally, the 3D fractal dimension (fd) of the simulated intralobar vascular tree was computed using an original implementation of classical box counting within the Mathematica scientific environment, as outlined in Supplementary Material under Computation of the Fractal Dimension of the Simulated Intralobar Vascular Tree section.

Validation procedures

Validation of the computational robustness of our IAN simulations with respect to the vascular anatomy of the human thyroid gland was performed by evaluating the 3D fd of the DLA structure with respect to an injection–corrosion cadaveric thyroid vascular cast (TVC). In brief, the gland of a human subject (male, age 59 years, BMI 24 kg/m²) obtained at autopsy from the Section of Forensic Medicine of the SBIBIT Department at the University of Parma (courtesy of Giovanni Lanzi), Italy, was injected through the superior thyroid artery (STA) and inferior thyroid artery (ITA) with a solution of 3% polyvinyl chloride covinyl acetate (Pevikon C 870; Sigma-Aldrich) diluted in acetone, to fill the IAN.

After resin polymerization, the thyroid tissue was digested by immersion in a solution of 5% HCL/pepsin for 48–72 h, leading to retrieval of the glandular IAN.^5,8,22,25 In a second step, the TVC was scanned using the microtomographic system Skyscan 1172 (Bruker mCT, Belgium) at the nominal resolution of 30 μm, through a rotation up to 180° with 0.4° rotation step and a frame averaging of 4, and reconstructed in three orthogonal planes using the Skyscan NRecon software (version 1.6.2.0, courtesy of Dr. Annapaola Parrilli, IOR, Bologna, Italy). Finally, the right lobe of the reconstructed TVC was selected, and a mean 3D fd was calculated using the Bruker CT-Analyser software (version 1.12.0) as previously applied.²³

An additional validation was achieved by comparing the percentage of lobar arteries belonging to each vascular order between computer reconstructions and literature data on the human thyroid gland.^30–32 In a first step, we calculated the absolute number of (a) node calibers (including roots) of 2nd order vessels, (b) node calibers of 3rd order vessels, and (c) node calibers of 4th order vessels up to an intermediate caliber of 0.025 mm (instead of the smallest caliber of 0.0125 mm; see the caliber intervals table in Supplementary Material).

This procedure was repeated for each simulation. The need for considering an intermediate caliber was due to excessive computational time; as a consequence, the number of 4th order nodes remained incomplete. Therefore, we estimated all of them (i.e., up to caliber 0.0125 mm) by doubling the number of computed 4th order nodes. This choice was based on the heuristic that a maximum of eight applications of the Roux–Murray Law (with exponent δ = 2.69; see formula (5) in Supplementary Material) were needed to span the entire 4th order interval of calibers, as opposed to the four applications expected to reach the caliber 0.025 mm.

Node calculations were achieved by performing an inverse scanning of the simulated tree, whose information was stored in four arrays organized in list form as follows:

- nodeList: containing the spatial coordinates of the nodes in the cluster,

- sizeList: containing the calibers of the nodes in the cluster,

- fatherList and sonList: containing integer indexes relating each node to its sons and father up to the 2nd order root.

Then, an indexing method was developed. Specifically, this is an algorithm with memory of performed steps where each node (or bifurcation) was tagged exploiting the list structure in Mathematica. This allowed for identification of familial relationships including ancestor and descendant particles, and their coordinates and calibers as well. In this way, information on nodes (including 2nd order roots) belonging to each vessel order was gathered for all orders.

In a second step, we developed a logical frame to estimate the number of 3rd and 4th order vessels in our simulations as follows:

(i)

The number of simulated 2nd order arteries was arbitrarily chosen at the minimum of 16 roots, corresponding to the minimal number of 2nd order arteries in a thyroid.^30–32

(ii)

Based on the number of computed 2nd order nodes (and roots), and having evidence that their values represented 16 real arteries (i.e., the same number arbitrarily chosen to start our simulations), we estimated that the mean number of algorithm iterations necessary to reconstruct 16 arteries per simulation was the ratio between the mean of 2nd order nodes (and roots) and the starting number of arterial roots (i.e., 16). This resulted in an average of 47 algorithm iterations per tree (equivalent to the number of 2nd order nodes/vessel) that spanned the entire caliber interval (of width equal to 0.34 mm) of 2nd order vessels.

(iii)

To obtain an estimate of the unknown number of 3rd order vessels in our simulations, we reasoned that if 47 algorithm iterations spanned the 2nd order interval width (0.34 mm), then the number of algorithm iterations to span the entire caliber interval (of width equal to 0.056 mm) of 3rd order vessels could be obtained as a proportion between these values, resulting in about eight algorithm iterations per tree on average (specifically, 7.76). Thus, the number of 3rd order estimated vessels resulted as the ratio between the mean simulated 3rd order nodes and these eight iterations.

(iv)

Finally, to obtain an estimate of the unknown number of 4th order vessels in our simulations, we reasoned that if an average of eight algorithm iterations spanned the entire 3rd order interval width (0.056 mm), then the number of algorithm iterations to span the entire caliber interval (of width equal to 0.016 mm) of 4th order vessels could be obtained as a proportion between these values, resulting in about two iterations per tree on average (specifically, 2.15). Thus, the number of 4th order estimated vessels resulted as the ratio between the mean estimated 4th order nodes and these two iterations.

In a third step, an estimate of the number of 2nd, 3rd, and 4th order vessels for an average human thyroid lobe was obtained taking into account the following anatomical data:^30–34

(a) The number of 1st order arteries collectively stemming from STA and ITA (four to five branches).

(b) The number of 2nd order vessels stemming from each 1st order artery (four to five branches).

(d) The anatomical evidence that, on average, each thyroid lobule is vascularized by a single 3rd order vessel.

(e) The number of 4th order vessels stemming from each 3rd order vessel (two to five branches).

Using these informations, minimal, mean, and maximal interval values for the number of vessels in each order were determined; each of these intervals was used to provide the percentages of vessels in a thyroid lobe. Finally, we calculated the distance (Δ%) between mean simulated and real vascular percentage values.

Statistical analysis

Differences in the 3D fd (either by order or averaged) of the simulated arterial trees were determined using one-way ANOVA. An r Pearson's correlation coefficient was computed for the overall mean 3D fd distribution of both simulated cases and vascular cast. Finally, differences in the distance (Δ%) between mean simulated and real vascular percentage values were calculated by initial arcsin transformation, and comparison of the transformed percentages using a z-value in a so-called simplest statistical test.³⁵ All results were considered statistically significant if p < 0.05.

Results

Figure 1 summarizes the steps for 3D reconstruction of the thyroid IAN in the three subjects studied. Lobar simulation started from a sagittal ultrasonographic image of the right thyroid lobe. After edge detection and spline approximation of its outlined contours, the lobe was reconstructed as a polygonal mesh, providing a closed space for growth of the IAN by the constrained DLA algorithm. Vascular simulations started from 16 vascular roots, each representing 1 of 16 real 2nd order thyroid arteries placed on the lobar surface keeping a reciprocal uniform distance.

FIG. 1.

Schematics of the steps for 3D reconstruction of the IAN in the three subjects studied (cases 1, 2, and 3). The procedure started from an ultrasonographic longitudinal section of their right thyroid lobe taken at its longest sagittal axis; (A) thyroid lobe profiles after edge detection and spline approximation, including their pseudosymmetry axis; (B) 3D reconstruction of the right lobe of each patient; (C–E) DLA growth of the IAN of the right thyroid lobe of each patient from 2nd to 4th order vessels, respectively. For mere rendering reasons, arteries stemming from roots located roughly above and below the barycentric plane parallel to the XY-plane were depicted with different colors, respectively. 3D, three-dimensional; DLA, diffusion-limited aggregation; IAN, intralobar arterial network. Color images are available online.

A progressive increase in the IAN density (from 2nd to 4th order vessels) occurred and, for mere rendering reasons, reconstructed arteries were depicted with different colors depending on roots located above and below the barycentric plane parallel to the XY-plane, respectively. Similar to normal anatomy, superior and inferior vessel groups tended to overlap at the level of the barycentric plane and, moving from ovoid-like (case 1) to conic-like (cases 2 and 3) lobar morphologies, a diffusion of superior vessels within the upper lobe pole occurred with respect to inferior vessels segregating in the lower lobe pole.

Figure 2 depicts the values of the 3D fd of the IAN in each of the simulated cases in comparison with that of a cadaveric case. A statistically significant difference emerged in the mean value of 3D fd for 2nd order vessels with respect to mean values of 3rd and 4th order vessels, as expected for a multifractal arterial pattern. However, a close degree of self-similarity in arterial geometry was found among the three IAN reconstructions for all order vessels, leading to absence of any statistically significant difference in their mean 3D fd values. These data linearly correlated in a statistically significant manner with the mean 3D fd value of the virtual microtomographic reconstruction of the cadaveric vascular cast.

FIG. 2.

(A) Values of the 3D fd of the IAN in each of the simulated cases analyzed. Note that a statistically significant difference occurs in the mean value of 3D fd for 2nd order vessels with respect to mean values of 3rd and 4th order vessels, suggesting a multifractal vascular pattern. In contrast, no statistically significant differences emerged between mean 3D fd based on all vascular orders, confirming consistency of the generative process of our algorithm. (B1) Human thyroid gland collected at autopsy. (B2) Corrosion cast of its IAN obtained after injection of synthetic resin. Note the presence of thyroid tissue in the upper pole of the right lobe, due to a still incomplete process of tissue digestion. (B3) Rendering of the reconstruction by microtomographic analysis of the resin-injected IAN of the right thyroid lobe, after completion of the tissue digestion process. (C) Linear regression analysis of the mean 3D fd of the IAN values of both simulated cases and vascular cast. Note the statistically significant correlation between all data. *p < 0.05 with respect to 3rd and 4th order vessels; NS, not statistically significant with respect to cases 2 and 3. fd, fractal dimension. Color images are available online.

Figure 3 summarizes results of the indexing method to calculate the number of roots and nodes of 2nd order arteries, and nodes of 3rd and 4th order arteries, in the simulations of the three euthyroid control cases. It also reports the interval of values for the number of 2nd to 4th order arteries in an average thyroid gland as deduced from literature data, and the differences in the proportion (percentage) of arteries at each vascular order between simulated and real cases.

FIG. 3.

(A) Number of roots and nodes of 2nd order arteries, and nodes of 3rd and 4th order arteries, as determined by the indexing method in the simulations of the IAN of our three euthyroid control cases. Values in the last column refer to the mean percentage of vessels at each artery order with respect to the total of the lobe arteries. (B) Estimated ranges of the number and percentages of arteries in the adult human thyroid lobe, as determined by anatomical literature data. Note the wide intervals, spanning from minimal to maximal through mean values (in bold). (C) Distribution of the difference in percentage between number of arteries in simulated versus real cases using mean artery values. The simulations did not yield statistically significant oscillations in the number of simulated vessels with respect to known anatomical data, suggesting that the generative principles of our algorithm robustly overlapped those of the native real vessels. min, inferior extreme of the interval; mean, average level of the interval; max, superior extreme of the interval; NS, not statistically significant.

Both in simulated and real cases, the number of arteries increased following a similar trend, with a vascular “explosion” (∼80 times) from 2nd to 3rd order, and a more limited ramification (∼12 times) from 3rd to 4th order. Differences in arterial percentages between simulated and real arteries were not statistically significant, confirming robustness of our algorithm in reproducing an IAN consistent with known anatomical data.

Figure 4 shows the flowchart of the procedure for simulation of the IAN of the human thyroid lobe. Phase 1 was set to describe the surface of the 3D thyroid lobe, giving rise to a 3D numerical model graphically rendered, whose volume could be estimated. The output of phase 1 provided a primary constraint for the 3D growth of the IAN, set in phase 2. The latter represented the core of the template, finalized to the generation of the IAN using a constrained DLA algorithm. The output of phase 2 constituted the set of data whose organization was exploited in phase 3, to compute the 3D fd of the IAN. This 3D fd worked as a control block to either proceed to the final output in STereoLithography or standard to layer (STL) format, or to break the flow and feedback to phase 2 for modification of DLA parameters.

FIG. 4.

Computational template summarizing the 3D modeling of the IAN of the adult human thyroid lobe. After phase 3, note the presence of two subsequent control blocks based on the 3D fd and indexing method. The latter measures the number of roots and nodes at each vascular order, to allow for estimation of number of vessels with respect to their anatomical reality. Both are decision switches for either feeding back into the DLA growth process or exiting to a computational format suitable for 3D bioprinting, like STL. Block of mean % of IAN branches and related arrows can be activated at request. STL, standard to layer. Color images are available online.

When 3D fd was considered acceptable, an additional control block could be activated (indexing) to count the number of simulated nodes (and roots) belonging to different vascular orders (optional phase 4), to provide a mean percentage of the IAN branches at different orders as an indicator of IAN model accuracy. This mean vascular percentage could be used to either proceed to the final STL output, or to break the flow and feedback to phase 2 to improve the DLA simulation.

Discussion and Conclusions

In this study, we developed an original computational template of numerical and rendering routines for simulating the 3D IAN of the human thyroid gland, based on a DLA algorithm. A peculiarity of our approach was to constrain the vessel growth within the 3D shape of a virtually reconstructed human thyroid lobe, with boundaries set through an homotopic-like parametrization involving splines and Bezier curves.^36,37

This choice (detailed in Supplementary Material) stemmed from previously tested building blocks developed to provide a robust and flexible framework yielding a reproducible structure to be eventually used for additive layer manufacturing/bioprinting of the stromal matrix or SVS of human organs, primarily the human thyroid gland.^5,7,21–27

Indeed, the 3D geometry of the SVS acts as a “photographic negative” of the macro-microscopic organ morphology, being per se capable of mechanochemically and topobiologically guiding cells during embryonic development and postnatal growth.^1,38 Therefore, we raised the possibility that the human thyroid lobe is organized as a topological “phase space” where growing thyrocytes may self-assemble driven by rules inherent in the 3D geometry of the SVS,^5–7 which acts as a fractal-like chaotic attractor.^7,9,21–25 Indeed, fractality is a property belonging to numerous human vascular networks.^39,40

Hence, we conceived a thyromorphic scaffold-bioreactor unit potentially scalable as a part of a microfluidic chip test system whose basic requirement resides in the strict similarity between the 3D geometry of the ex situ engineered/bioprinted SVS and the natural thyroid lobe SVS.^5,8,9 We recently substantiated this perspective showing that a decellularized organomorphic, rodent, and pig SVS is suitable for in vitro reconstruction of lung, liver, thymus, adrenal, anterior pituitary, cerebral cortex, pancreas, kidney, skin, and bone tissues by recellularization.^8,41–43

Consistently, we simulated a patient-tailored SVS constrained within the geometric boundaries of the virtually reconstructed lobar “phase space,” and modeled as an hybrid fractal network with a “deterministic boost” provided by average values of some standard vascular parameters (for details see Supplementary Material). Exploiting the evidence that in the normal thyroid gland there is spatial uniformity of the lobe fitting an ovoid-like or conic-like geometry,^5,6,9 we reproduced the 3D lobar shape from a single ultrasonographic lobar image (for details see Supplementary Material). We also calculated the volume of these 3D lobar reconstructions using either our algorithm or, for comparison, different and more routinely used calculations.^44–48 All volumes resulted substantially coincident (data not shown).

Once the DLA was iterated up to 4th order arteries, we calculated the 3D fd (i.e. fractal dimension) values of this arterial distribution using an ad hoc box counting method (see Supplementary Material). Fractal dimension values resulted statistically similar in all three reconstructed thyroid lobes, and linearly correlated with those of a vascular cast from a cadaveric thyroid used as an internal control. However, comparison between mean 3D fd of 2nd with respect to 3rd and 4th order arteries depicted a statistically significant difference, suggesting multifractality in the branching pattern, as occurring in a number of vascular arborizations of human viscera.⁴⁹

To substantiate the similarity between DLA simulations and current anatomical data on thyroid lobe vessels,^30–33 we developed a heuristic based on a control routine counting simulated roots and nodes by vessel order (so-called indexing method) coupled to an estimate of the number of algorithmic iterations needed to “utilize” the entire interval of real vessel calibers (as deduced from known literature data) in each simulated SVS. In this manner, an estimate of the number of 2nd to 4th order simulated vessels was reached (for details on this procedure see Supplementary Material).

As a result, we obtained a percentage of 2nd to 4th order vessels not statistically different from that predictable by current anatomical data (details on this frame are in Supplementary Material). However, beyond automation of the indexing method, a patient-tailored SVS simulation by DLA might benefit of a machine learning procedure aimed at classifying and recognizing vascular patterns corresponding to those acquired by in vivo imaging techniques on each single subject. Similarly, DLA has recently been combined with neural network algorithms to improve diagnosis of brain injury of individual patients.⁵⁰ Whether this optimization should proceed by supervised inputs or reinforcement mechanisms⁵¹ remains to be tested.

Finally, we developed an operational flowchart of the entire process aimed at the 3D virtual reconstruction of the thyroid SVS, savable in STL format for its prototyping with biomaterials.^{9,21,22,24–26} Use of an STL-based 3D blueprint of ramified vascular structures is a well-known requirement for their fabrication by extrusion bioprinting^52,53 and material jetting technologies.⁵⁴ Consistently, in previous studies based on preliminary blocks of our computational template, we tested the suitability of a multi jet modeling approach to fabricate with Ultraviolet (UV)-curable resins a single branch prototype of the human thyroid SVS.

Results showed that minimal calibers of printed endings varied from 0.11 to 0.41 mm, with an excess dimensional deviation of 0.01–0.03 mm with respect to the STL simulation,²² thus depicting acceptable limits of dimensional accuracy for 2nd order vessels. Recent development of material extrusion and jetting technologies, like layer-by-layer ultraviolet-assisted extrusion-based and drop-on-demand inkjet bioprinting with hydrogels,^54,55 promises to reach dimensional limits and accuracy compatible with those requested by 3rd and 4th order SVS branches. Eventually, when in vivo tomographic technologies will allow for collection of volumetric images of the intrinsic vascularity of human organs, volumetric additive manufacturing by tomographic images polymerization will prove ideal for patient-tailored reconstructions,^56–58 including thyroid SVS.

To model fluid dynamics in mammalian vessels, a number of computational approaches relied upon fractal geometry^59–61 coupled to vascular permeability, pressure drop and difference, Poiseuille's equation, and Reynold's number.⁶² However, our primary focus was on a geometric structural similarity with the anatomical reality of the native thyroid SVS, at a scale level above that of capillary beds requiring different network laws.^39,63,64 Since embryonic and fetal development of vascular systems leads to an adult vascular geometry through a hemodynamic-induced remodeling action,^65,66 it may be assumed that the 3D geometry of an adult SVS reasonably retains structural features compatible with the organ-specific fluid dynamics.⁹

For these reasons, we predicted that the low/medium geometric constraints offered by our 2nd to 4th order IAN would allow for perfusion dynamics of culture media, favoring maintenance and self-aggregation of metastable toroidal patterns of seeded cells into luminous modules (follicle-like structures). Later, cocultured endothelial precursors would be expected to give rise to a dense interwoven perifollicular capillary bed (so-called thyroid baskets) whose highly geometric constraints could lead to self-limited follicular inflation, and arrest of glandular growth within the limits of the thyroid lobe boundaries.^5,9

In conclusion, we believe that the computational template introduced here might prove efficient to engineer an organomorphic scaffold reproducing the 3D macro-microscopic architecture of the adult human thyroid lobe. This would offer a powerful tool for its ex situ reconstruction using human adult thyroid stem cells and endothelial precursors, up to fabrication of a thyromorphic microfluidic chip test system. In addition, relying on standard morphometric variables such as vessel order and size, bifurcation law, number of branches per vessel order, and fd, it provides a general concept for SVS fabrication, suitable to ex situ bioprinting of various human organs upon clinical request.

Footnotes

Acknowledgments

Part of this study was developed under the tenure of the 2020–2023 DIMEC UNIPR-CMG Collaboration Agreement signed on October 8, 2020, as a part of the research program: Outpatient Diagnostic—Therapeutic Quality and National and International Guidelines on Endocrine-Metabolic Disorders. The authors are also grateful to the Medical Center Galliera in San Venanzio di Galliera, Bologna, Italy, for providing clinical and technical resources to pursue ultrasonographic studies of the thyroid gland in human subjects.

Authors' Contributions

G.S. and M.S. contributed to conceptualization and analysis of numerical mathematical procedures, ideation and structure of the computational kernel, and software resources; F.B., G.C., and S.M. were involved in selection and anatomical recognition of ultrasonographic thyroid images, collection of thyroid cadaveric anatomical specimens, injection–corrosion procedure for vascular identification, and application of routines for statistical calculations; R.T. carried out research conceptualization, development and supervision of the entire scientific program, collection of clinical ultrasonographic images, evaluation of appropriateness of all methodologies (imaging, anatomical, mathematical, geometric, and statistical), data evaluation and heuristic of estimates, article writing, editing, and supervision including English language.

Disclosure Statement

No competing financial interests exist.

Funding Information

Part of this research was conducted under the support of the UNIPR-FIL 2020–2022 Fund and Emilia-Romagna Sisma Ripopolamento 2/2019 Fund at the Medical Center Galliera in San Venanzio di Galliera, Bologna, Italy. Ensuing technologies contributed to improving the studies under the tenure of the Horizon 2020 SCREENED research program, no. 825745 ().

Supplementary Material

References

Edelman

GM.

Topobiology: An Introduction to Molecular Embryology. Basic Book, Inc.: New York, USA; 1988.

Chuong

, Wu

, Plikus

, et al. Engineering stem cells into organs: Topobiological transformations demonstrated by beak, feather, and other ectodermal organ morphogenesis. Curr Top Dev Biol, 2006; 72:237–274.

Fagman

, Andersson

, Nilsson

The developing mouse thyroid: Embryonic vessel contacts and parenchymal growth pattern during specification, budding, migration, and lobulation. Dev Dyn, 2006; 235:444–455.

Fagman

, Nilsson

Morphogenetics of early thyroid development. J Mol Endocrinol, 2011; 46:R33–R42.

Toni

, Della Casa

, Spaletta

, et al. The bioartificial thyroid: A biotechnological perspective in endocrine organ engineering for transplantation replacement. Acta Biomed, 2007; 78 (Supp 1):129–155.

Toni

, Della Casa

, Bodria

, et al. A study on the relationship between intraglandular arterial distribution and thyroid lobe shape: Implications for biotechnology of a bioartificial thyroid. Ann Anat, 2008; 190:432–441.

Toni

, Tampieri

, Zini

, et al. Ex situ bioengineering of bioartificial endocrine glands: A new frontier in regenerative medicine of soft tissue organs. Ann Anat, 2011; 193:381–394.

Sprio

, Sandri

, Iafisco

, et al. Biomimetic Biomaterials in Regenerative Medicine. In: Biomimetic Biomaterials, Structure and Applications. (Ruys AJ. ed.) Whoodhead Publishing: Cambridge, UK, 2013; pp. 3–45.

Toni

, Bassi

, Barbaro

, et al. Bioartificial endocrine glands: At the Cutting Edge of Translational Research in Endocrinology. In: Bio-Inspired Regenerative Medicine: Material, Processes, and Clinical Applications. (Sprio S, Tampieri A. eds.) Pan Standford Publishing: Singapore, 2016; pp. 357–387.

10.

Strusi

, Zini

, Dallatana

, et al. Ex situ bioengineering of the rat thyroid using as a scaffold the three dimensional (3D) decellularized matrix of the glandular lobe: Clues to the organomorphic principle. Ital J Anat Embryol, 2011; 116 (Suppl 1):S180.

11.

Strusi

, Zini

, Dallatana

, et al. Endocrine bioengineering: Reconstruction of a bioartificial thyroid lobe using its three-dimensional (3D) stromal/vascular matrix as a scaffold. Endocr Abstr (94th Meeting, Huston), 2012; 29:1586.

12.

Toni

, Strusi

, Zini

, et al. Bioengineering of the thyroid lobe: Use of its stromal/vascular matrix as a scaffold for ex situ reconstruction. Endocr Rev, 2012; 33 (03_Meeting Abstract):OR26–3.

13.

Alfieri

, Barbaro

, Consolini

, et al. A targeted mass spectrometry method to screen collagen types I–V in the decellularized 3D extracellular matrix of the adult male rat thyroid. Talanta, 2019; 193:1–8.

14.

Ogundipe

VML

, Plukker

JTM

, Links

, et al. Thyroid gland organoids: Current models and insights for application in tissue engineering. Tissue Eng Part A, 2022; 11–12:500–510.

15.

Bulanova

, Koudan

, Degosserie

, et al. Bioprinting of a functional vascularized mouse thyroid gland construct. Biofabrication, 2017; 9:034105.

16.

Kurmann

, Serra

, Hawkins

, et al. Regeneration of thyroid function by transplantation of differentiated pluripotent stem cells. Cell Stem Cell, 2015; 17:527–542.

17.

Ogundipe

VML

, Groen

, Hosper

, et al. Generation and differentiation of adult tissue-derived human thyroid organoids. Stem Cell Reports, 2021; 16:913–925.

18.

Pan

, Li

, Fang

, et al. Regeneration of a bioengineered thyroid using decellularized thyroid matrix. Thyroid, 2019; 29:142–152.

19.

Weng

, Chen

, Xie

, et al. Rabbit thyroid extracellular matrix as a 3D bioscaffold for thyroid bioengineering: A preliminary in vitro study. Biomed Eng Online, 2021; 20:18.

20.

Aleshcheva

, Bauer

, Hemmersbach

, et al. Scaffold-free tissue formation under real and simulated microgravity condition. Basic Clin Pharmacol Toxicol, 2016; 119:26–33.

21.

Bassoli

, Denti

, Gatto

, et al. New Approaches to Prototype 3D Vascular-Like Structures by Additive Layer Manufacturing. In: Innovative Developments in Virtual and Physical Prototyping. (Bartolo PJ, Soares De Lemos AC, Oliveria Tojeira AP, et al. eds.) CRC Press, Taylor & Francis: London, UK, 2012; pp. 35–42.

22.

Bassoli

, Denti

, Gatto

, et al. A combined additive layer manufacturing (ALM)/indirect replication method to prototype 3D vascular-like structures of soft tissue and endocrine organs. Virtual Phys Prototyp, 2012; 7:3–11.

23.

Parrilli

, Zini

, Spaletta

, et al. Micro-CT of 3D Vascular Structure: Clues for Innovative Scaffolds in Organ Engineering. In: Proceeding Bruker MicroCT User Meeting. (Brucker Instruments. ed.) Hasselt, Belgium, 2013; pp. 135–139.

24.

Toni

, Zini

, Barbaro

, et al. Bioartificial Endocrine Organs: A Translational Perspective in Regenerative Endocrinology. In: Proceedings of the MiMe Materials in Medicine (CNR–ISTEC. eds). CNR–ISTEC: Faenza, Italy, 2013; pp. 102-MC-O.

25.

Bassoli

, Denti

, Gatto

, et al. A Planar Fractal Analysis of the Arterial Tree of the Human Thyroid Gland: Implications for Additive Manufacturing of 3D Ramified Scaffolds. In: High Value Manufacturing. Advanced Research in Virtual and Rapid Prototyping. (Bartolo PJ, Soares De Lemos AC, Pereira AMH, et al. eds.) CRC Press, Taylor & Francis: London, 2014; pp. 423–429.

26.

Bassoli

, Denti

, Gatto

, et al. Towards additive manufacturing of ramified scaffolds of the thyroid vascular system: A preliminary fractal analysis. Int J Mech Eng Technol, 2018; 9:429–437.

27.

Spaletta

Reconstruction in space and visualization of a planar image: A mathematical and computational introduction. Acta Biomed, 2007; 78(Supp 1):26–31.

28.

Wolfram

The Mathematica Book, 5th Edition. Wolfram Media: Champaign, IL, USA, 2003.

29.

Witten

Jr , Sander

. Diffusion-limited aggregation, a kinetic critical phenomenon. Phys Rev Lett, 1981; 47:1400–1403.

30.

Major

RH.

Studies on the vascular system of the thyroid gland. Am J Anat, 1909; 9:475–492.

31.

Wangesteen

OH.

The blood supply of the thyroid gland with special reference to the vascular system of the cretin goiter. Surg Gynecol Obstet, 1929; 48:613–628.

32.

Papadatos

Some anatomo-radiological observations concerning the changes in thyroid arteries which occur with senility. Anat Anz, 1981; 150:212–225.

33.

Murray

The Thyroid Gland. In: Functional Endocrine Pathology, Vol. 1. (Kovacs K, Asa S. eds.) Blackwell: Oxford, 1991; pp. 293–374.

34.

Brown

, AL-Moussa

, Beck

. Histometry of normal thyroid in man. J Clin Pathol, 1986; 39:475–482.

35.

Pocock

SJ.

The simplest statistical test: How to check for a difference between treatments. BMJ, 2006; 332:1256–1258.

36.

Whitehead

GW.

Elements of Homotopy Theory. In: Homotopy, Encyclopedia of Mathematics. (Hazewinkel M. ed.) Springer: Berlin, 1978; p. 228.

37.

Schoenberg

JT.

Contributions to the problem of approximation of equidistant data by analytic functions. Part A: On the problem of smoothing or graduation. A first class of analytic approximation formulae; Part B: On the problem of osculatory interpolation. A second class of analytic approximation formulae. Quart Appl Math, 1946; 4:45–99 and 112–141.

38.

Ingber

DE.

Mechanical control of tissue growth: Function follows form. Proc Natl Acad Sci U S A, 2005; 102:11571–11572.

39.

Zamir

On fractal properties of arterial trees. J Theor Biol, 1999; 197:517–526.

40.

Zamir

Arterial branching within the confines of fractal L-system formalism. J Gen Physiol, 2001; 118:267–275.

41.

Toni

, Di Conza

, Barbaro

, et al. Microtopography of immune cells in osteoporosis and bone lesions by endocrine disruptors. Front Immunol, 2020; 11:1737.

42.

Asnaghi

, Barthlott

, Gullotta

, et al. Thymus extracellular matrix-derived scaffolds support graft-resident thymopoiesis and long-term in vitro culture of adult thymic epithelial cells. Adv Funct Mater, 2021; 2010747:1–15.

43.

Remaggi

, Barbaro

, Di Conza

, et al. Decellularization detergents as methodological variables in mass spectrometry of stromal matrixes. Tissue Eng Part C, 2022; 28:148–157.

44.

Hegedus

, Perrild

, Poulsen

, et al. The determination of thyroid volume by ultrasound and its relationship to body, weight, age, and sex in normal subjects. J Clin Endocrinol Metab, 1983; 56:260–263.

45.

Szebeni

, Beleznay

New simple method for thyroid volume determination by ultrasonography. J Clin Ultrasound, 1992; 20:329–337.

46.

Miccoli

, Minuto

, Orlandini

, et al. Ultrasonography estimated thyroid volume: A prospective study about its reliability. Thyroid, 2006; 16:37–39.

47.

Ying

, Yung

, Ho

Two dimensional ultrasound measurement of thyroid gland volume: A new equation with higher correlation with 3D ultrasound measurement. Ultrasound Med Biol, 2007; 34:56–63.

48.

Slapa

, Jakubowski

, Slowinska-Srzednicka

, et al. Advantages and disadvantages of 3D ultrasound of thyroid nodules including thin slice volume rendering. Thyroid Res, 2011; 4:1–12.

49.

Zamir

Fractal dimensions and multifractality in vascular branching. J Theor Biol, 2001; 212:183–190.

50.

Karaca

, Cattani

, Karabudak

ANN Classification of MS Subgroups with Diffusion Limited Aggregation. In: Computational Science and Its Applications—ICCSA 2018 (Gervasi O, Murgante B, Misra S, et al. eds.), LNTCS 10961. Springer International Publishing AG: Cham, CH, 2018; pp. 121–136.

51.

Kim

, McKee

, Fontenot

, et al. Engineering tissue fabrication with machine intelligence: Generating a blueprint for regeneration. Front Bioeng Biotechnol, 2020; 7:443.

52.

Jiang

, Munguia-Lopez

, Flores-Torres

, et al. Extrusion bioprinting of soft materials: An emerging technique for biological model fabrication. Appl Phys Rev, 2019; 6:011310.

53.

Zhang

, Haghiashtiani

, Hübscher

, et al. 3D extrusion bioprinting. Nat Rev Methods Primers, 2021; 1:75.

54.

, Liu

, Pei

, et al. Inkjet bioprinting of biomaterials. Chem Rev, 2020; 120:10793–10833.

55.

Zhuang

, Long Ng

, An

, et al. Layer-by-layer ultraviolet assisted extrusion-based (UAE) bioprinting of hydrogel constructs with high aspect ratio for soft tissue engineering applications. PLoS One, 2019; 14(6):e0216776.

56.

Kelly

, Bhattacharya

, Heidari

, et al. Volumetric additive manufacturing via tomographic reconstruction. Science, 2019; 363:1075–1079.

57.

Long Ng

, Min Lee

, Zhou

, et al. Vat polymerization-based bioprinting: Process, materials, applications and regulatory challenges. Biofabrication, 2020; 12:022001.

58.

, Mille

, Robledo

, et al. Recent advances in formulating and processing biomaterial inks for vat polymerization-based 3D printing. Adv Healthc Mater, 2020; 9(15):e2000156.

59.

Rybaczuk

, Kediza

, Blaszczyk

Fractal description of cerebellum surface during fetal period. Folia Morphol, 1996; 55:434–436.

60.

Rybaczuk

, Kediza

, Andrzejak

Fractal dimension of human brain cortex vessel during the fetal period. Med Sci Monit, 2002; 8:46–51.

61.

Rybaczuk

, Kediza

, Paradowsky

Fractal characteristics of brain vessel microangioarchitecture during the fetal period. Med Sci Monit, 2002; 8:145–152.

62.

Jayalalitha

, Shanthoshini Deviha

, Uthayakumar

Fractal model of blood flow in cardiovascular system. Computers Biol Med, 2008; 38:684–693.

63.

Sherman

On connecting large vessels to small. The meaning of Murrays law. J Gen Physiol, 1981; 78:431–453.

64.

Lorthois

, Cassot

Fractal analysis of vascular networks: Insights from morphogenesis. J Theor Biol, 2010; 262:614–633.

65.

Risau

Mechanisms of angiogenesis. Nature, 1997; 386:671–674.

66.

le Noble

, Moyon

, Pardanaud

, et al. Flow regulates arterial-venous differentiation in the chick embryo yolk sac. Development, 2004; 131:361–373.

67.

Canuto

, Hussaini

, Quarteroni

, et al. Spectral Methods: Fundamentals in Single Domains. Springer Verlag: Berlin and Heildelberg, Germany, 2006.

68.

Kurtz

, Sandau

, Christ

On the bifurcation of blood vessels—Wilhelm Roux's Doctoral Thesis (Jena 1878)—A seminal work for biophysical modelling in developmental biology. Ann Anat, 1997; 179:33–36.

69.

Liu

, Chen

, Wang

, et al. An improved differential box-counting method to estimate fractal dimensions of gray-level images. J Vis Commun Image Represent, 2014; 25:1102–1111.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.59 MB

0.00 MB