The fractal spatial distribution of pancreatic islets in three dimensions: a self-avoiding growth model

Improving glucose homeostasis by regenerating the insulin-secretion capability in diabetics is a major medical priority. Insulin is secreted by beta cells in the endocrine pancreas. A normal beta-cell function depends on the collective behavior of endocrine cells organized in the islets of Langerhans. Thus, apart from the crucial question of increasing beta-cell mass, the process of formation and function of the islet distribution in a normal pancreas is of topical interest in the light of rapid advances in bioengineering.

The fractal geometry of self-similar structures is observed in many organs such as the lung, kidney, brain and vasculature [1–5]. Such structures may originate from developmental processes using basic building rules repeatedly, as confirmed in lung formation [6]. A defining property of fractals is the scale independence of self-similar patterns. The space filling property is measured by the fractal dimension, a quantity that can take non-integer values. The methods for measuring fractal dimensions from biological images have been reviewed [7]. The fractal dimensions, D₂, measured usually from two-dimensional (2D) organ images are 1.8 in the bronchial tree [1], 1.6 in the renal arterial tree [2], 1.5 in the cerebral cortex [3], 1.3 in the pial vasculature [4] and 1.7 in the retinal vasculature [5]. In addition to organ structures, the spatial distributions of islets on pancreas sections have been reported to follow a fractal geometry with D₂ = 1.5 in guinea pigs [8] as well as in dogs, pigs, goats and cows [9].

These fractal dimensions were computed from 2D images of organs. Recent advances in imaging methods, however, have provided three-dimensional (3D) information, e.g., the 3D fractal dimension of the lung, D₃ = 2.3–2.5 [10]. In particular, optical projection tomography (OPT) has provided 3D spatial distributions of islets in intact mouse pancreata [11]. In addition to such advances in imaging methods, the availability of computational resources has enabled practical advances in automated image analysis [12]. Applying these to human tissues has allowed the examination of the geometry of islet distributions in the entire human pancreata, although the data obtained from large human organs are still limited to 2D sections [13]. In this study, we examine the spatial organization of islets in 3D mouse pancreata and 2D human pancreatic sections using data obtained with these advances.

It has been speculated that the fractal geometry of the spatial islet distribution is associated with vascularization in the pancreas [8]. Here we explicitly examine this hypothesis with a self-avoiding growth model and show that the spatial fractal distribution of islets may result from self-avoidance in the branching process of vascularization.

The model starts with a single node (islet) at a 3D position, y₁. At the next growth step, the node can generate (branch into) a new node at y₂ = y₁ + Δy, where Δy is a random 3D vector with |Δy| < l. Here l is a distance scale obtainable from the mean distance between the nearest islets in experiments. In the following step, every node including the mother and daughter nodes can generate new nodes. In the growth process, we introduce self-avoidance to inhibit close overlaps between nodes. A standard mechanism for the directionality of growth processes is the diffusion of special molecules such as chemo-attractants and chemo-repellents, e.g., in axon guidance [14]. Using the idea of chemo-repulsion, we assume that every node at y_i produces a chemo-repellent factor, which diffuses out from the existing nodes and finally degrades. Given N nodes, the potential (concentration) of the chemo-repellents, Φ(x, t), at position x at time t can be described by

$$\begin{equation} \partial _t \Phi = \sum _{i=1}^N \delta ({\bf x}-{\bf y}_i) + D \partial _{\bf x}^2\Phi -d \Phi , \end{equation} \tag{ 1 }$$

where D and d are diffusion and degradation constants. Since the diffusion and degradation are expected to occur much faster than the growth process, we consider a stationary state of the above equation and obtain a spatial gradient of the chemo-repellents. We can solve this differential equation by the Fourier transform:

$$\begin{equation} \Phi ({\bf k}) = \frac{1}{D{\bf k}^2 + d} \sum _{i=1}^N {\rm e}^{-{\unicode{x131}} {\bf k} \cdot {\bf y}_i}. \end{equation} \tag{ 2 }$$

Finally, the inverse Fourier transform of Φ(k) gives

$$\begin{equation} \Phi ({\bf x}) = \sum _{i=1}^N \frac{\lambda }{|{\bf x}-{\bf y}_i|}\, {\rm e}^{-|{\bf x}-{\bf y}_i|/\lambda }, \end{equation} \tag{ 3 }$$

where $\lambda =\sqrt{D/d}$ is an effective distance of the chemo-repulsion from the existing nodes. In this study, we set λ = l to match the effective distance as the mean minimal distance between islets. It is of interest that Φ(x) has the form of the Yukawa potential in particle physics. Therefore, a new node can be generated in a unit volume around a point, x, near the existing ith node (|x − y_i| < λ) with the probability density

$$\begin{equation} P({\bf x}) = \mathcal {A} \Phi ^{-m} ({\bf x}), \end{equation} \tag{ 4 }$$

which is assumed to be proportional to the inverse of the chemo-repellent gradient from the existing nodes. The parameter m determines how strongly the chemo-repulsion affects the growth process with a corresponding normalization constant, $\mathcal {A}$. When no additional nodes exist near the ith node, the normalization constant is approximated as $\mathcal {A}^{-1}(m) \approx \int _0^1 r^m {\rm e}^{mr} {\rm d}r$. However, depending on the exponent m, a few nodes exist together within the effective distance λ from the ith node. For example, there are ten and five nodes on average within the effective distance for m = 0.5 and 3, respectively. In practice, it is not necessary to compute the exact normalization constant $\mathcal {A}$, which requires a heavy computation at every addition of new nodes. Instead, we use a practical normalization constant, $\mathcal {A}^{\prime }( {<} \mathcal {A}$), that equally reduces P(x) at every position, x. We have checked that the growth process is independent of the specific value of $\mathcal {A}^{\prime }$ if it is small enough (e.g. $\mathcal {A}^{\prime } = 0.01 \mathcal {A}$), while it must also be large enough not to reject every trial of node addition.

The practical algorithm for adding a new node to the existing N nodes is as follows.

• 1.  
  Calculate the chemo-repellent gradient Φ(x) in equation (3), generated by the existing N nodes.
• 2.  
  Randomly select one node (e.g. the ith node) among the N nodes.
• 3.  
  Randomly pick a coordinate near the ith node, x = y_i + δy, where y_i is the coordinate of the ith node and δy is a random 3D vector with |δy| < l.
• 4.  
  Accept the addition of the new node with a probability, P(x), in equation (4).
• 5.  
  Iterate these steps to add more nodes.

3.1. Fractal spatial distribution of pancreatic islets

OPT imaging captured the 3D spatial distribution of the islets scattered in the exocrine tissue of the mouse pancreas (figure 1(A)). Based on the centroid coordinate of every islet, we calculated the average number 〈n〉 of neighboring islets within a distance r. This allowed us to determine the correlation dimension D [15], a type of fractal dimension, from the relation 〈n〉∝r^D. For example, when points are homogeneously distributed in a line, surface and volume, their fractal dimensions become 1, 2 and 3, respectively, as expected in the Euclidean space. The logarithmic plot of 〈n〉 versus r showed the expected power-law behavior in 3D (figure 1(B)), first observed in the 2D pancreatic sections in guinea pigs with the fractal dimension D_2S = 1.5 [8]. Independent of pancreatic lobes (splenic, duodenal and gastric lobes), the 3D fractal dimension was D₃ = 2.1 (table 1). The box-counting method gave a similar 3D fractal dimension, 2.03 ± 0.24 (P = 0.3).

Zoom In Zoom Out Reset image size

Table 1. Fractal dimensions of the spatial islet distributions in mouse and human pancreata. The fractal dimensions of islet distributions in three-dimensional volume (D₃) and two-dimensional projections (D_2P) and sections (D_2S) of the pancreas were calculated based on islet centroid coordinates. N_total represents the total number of islets in the pancreas volume, while N_section represents the number of islets on the pancreas section. Mean ± s.d. (n samples).

Species	Age	n	Region	Ntotal	Nsection	D3	D2P	D2S
			Splenic	2036 ± 269	102 ± 18	2.06 ± 0.06	1.60 ± 0.03	1.28 ± 0.04
Mousea	8 weeks	5	Duodenal	1672 ± 198	100 ± 13	2.09 ± 0.09	1.60 ± 0.05	1.27 ± 0.02
(C57BL/6)			Gastric	778 ± 73	78 ± 17	2.13 ± 0.09	1.62 ± 0.03	1.37 ± 0.08
	All	15	–	–	–	2.10 ± 0.08	1.61 ± 0.04d	1.31 ± 0.07e
Mouseb	1 day	8	Whole	1204 ± 250	–	–	1.57 ± 0.04	-
(MIP-GFP)	3 weeks	3	Whole	3335 ± 709	–	–	1.54 ± 0.04	-
	All	11	–	–	–	–	1.56 ± 0.04d	-
Humanc	50 ± 20 years	14	Whole	–	370 ± 156	–	–	1.53 ± 0.07e

^aPancreata from female mice were isolated, stained for insulin and subjected to optical projection tomography as described in our previous study [24]. The 3D islet centroid coordinates were extracted using the Imaris software. For the estimation of the 2D projection and section fractal dimensions, the projection plane was selected by rotating the pancreas in 3D to maximize the averaging distance between islets on the projected plane, and the section plane along with the given projection plane was selected to include the maximal number of islets within a section depth of 100 μm. ^bPancreata were excised from MIP-GFP mice in which beta cells are genetically tagged with green fluorescent protein under the control of the mouse insulin 1 promoter [36]. Fixed and cleared pancreata were then placed between a slide glass and cover slip. To capture entire islets including those that are not in perfect focus, an epifluorescent configuration with the uncertainty of Z depth was used with a confocal microscope. The islet centroid coordinates were extracted using the ImageJ software. Detailed methods are described in our previous study [37]. ^cAutomated image analysis was used to examine centroid coordinates of islets on pancreatic sections from human autopsies [13]. ^dThe fractal dimensions of projection D_2P showed no significant difference between C57BL/6 and MIP-GFP mice. ^eThe fractal dimensions of sections D_2S showed a significant difference between C57BL/6 mice and humans (P < 0.01)

The spatial islet distribution information in the pancreas is usually extracted from 2D pancreatic sections, especially for large animals including humans. To examine the method dependence, we compared islet distributions in the 3D, 2D projections and 2D sections of the full 3D data (figure 2(A)). Note that random sections of the non-symmetric 3D structure in figure 1(A) could result in large variations in the fractal dimension measurement. Therefore, we consistently chose the unique 2D projection plane and section that maximize the average distances between islets and the number of islets within the section, respectively, by rotating the 3D pancreas. Although all three methods gave power-law behaviors in the spatial distribution, they had different fractal dimensions (figure 2(B) and table 1). The 3D fractal dimension has been extrapolated as D₃ = D_2S + 1 by adding a unit co-dimension based on the fractal dimension D_2S measured on the 2D pancreatic section [8]. Our direct observation demonstrated that D₃ was different from D_2S + 1 (P < 0.01), suggesting that islet distributions are not symmetric in space. This was obvious in the 3D OPT image (figure 1(A)).

Zoom In Zoom Out Reset image size

We also examined islet distributions in intact mouse pancreas with uncertainty in Z depth, which corresponds to the 2D projection. The fractal dimension did not change with age (table 1), which is consistent with the age independence found in guinea pigs [9].

Furthermore, the recent availability of human tissues allowed us to examine the spatial islet distribution in human autopsies [13]. We analyzed the fractal dimension with pancreatic sections from adult pancreata. They also followed the power law with the fractal dimension D_2S = 1.53 ± 0.07 (table 1), which is very close to 1.5, as reported for other smaller mammals [9]. Note that a relatively smaller animal, mouse, showed a significantly smaller fractal dimension, D_2S = 1.3 (table 1). In our model, this corresponds to a smaller value of m.

3.2. The emergence of fractal geometry

Pancreatic islets are co-localized with pancreatic ducts [16] and blood vessels [17] in the pancreas. Thus, it has been speculated that the fractal geometry of spatial islet distributions is associated with the vascular structure of pancreatic ducts [8, 18]. In pancreatic development, the primitive endodermal epithelium differentiates into the pancreatic duct and the endocrine and exocrine cells [19]. Note that the clusters of endocrine cells finally form islets. Blood vessels play a critical role for organ development by providing not only oxygen and nutrients, but also inductive signals for endocrine cell differentiation [20]. On the other hand, endocrine cells produce the vascular endothelial growth factor, which attracts endothelial cells of blood vessels. Therefore, the organogenesis of endocrine pancreas results from the mutual interaction between the endocrine cells and blood vessels. Surprisingly, however, Magenheim et al have recently reported that blood vessels reduce the branching and differentiation of epithelial cells [21]. The dual effect of inducing and inhibiting endocrine cell differentiation by blood vessels may lead to a self-avoiding formation of islets in the pancreas.

Here, we developed an organ growth model considering self-avoidance. The model focused on islet formation without introducing the epithelium and blood vessels explicitly. Thus islets, represented as nodes in space, were assumed to produce a chemo-repellent potential, Φ, to implicitly incorporate the self-avoiding effect resulting from the mutual interaction between the islets and blood vessels. In the model, new islets (nodes) were generated from the existing nodes with a probability P∝Φ^−m. The parameter m controls the degree of self-avoidance. For example, a vanishing m amounts to neglecting self-avoidance, while a large m limits the growth to locations furthest from the existing nodes. Figure 3 shows examples of the self-avoiding growth with various values of m. With a negligible self-avoidance (m = 0.5; figure 3(A)), nodes were densely distributed with a high fractal dimension (D₃ = 2.74 ± 0.04). On the other hand, with a relatively strong self-avoidance (m = 3; figure 3(C)), nodes were sparsely distributed with a small fractal dimension (D₃ = 1.88 ± 0.05). Finally, the self-avoiding growth model with m = 1.5 could produce a node distribution pattern resembling the islet distribution in the pancreas (figure 3(B)). The fractal dimensions of the 3D, 2D projections and 2D sections of the model were consistent with the fractal dimensions of islet distributions, when the node number is truncated to 2000 as the total islet number in adult mice (tables 1 and 2). The self-avoiding growth model could generate node distributions with distinct fractal dimensions (figure 3(D) and table 2). Note that during the growth process with a fixed m, the fractal dimension generally increased until the node number became large enough to make finite size effects negligible (table 2).

Zoom In Zoom Out Reset image size

Table 2. Fractal dimensions in the self-avoiding growth model. The degree of self-avoidance produced different fractal structures. A larger m produces stronger self-avoidance. The model had finite size effects on the fractal dimension. As the model with m = 1.5 grows from 1000 to 5000 nodes, its fractal dimensions of the 3D, 2D projections and sections approach equilibrium. Mean ± s.d. (n = 10).

m	Ntotal	Nsection	D3	D2P	D2S
0.5	2000	134 ± 17	2.47 ± 0.04	1.73 ± 0.02	1.45 ± 0.10
1.0	2000	110 ± 8	2.21 ± 0.05	1.64 ± 0.04	1.29 ± 0.09
1.5	2000	87 ± 12	2.06 ± 0.04	1.60 ± 0.07	1.22 ± 0.08
2.0	2000	84 ± 12	2.00 ± 0.05	1.59 ± 0.07	1.14 ± 0.08
3.0	2000	86 ± 11	1.88 ± 0.05	1.57 ± 0.05	1.11 ± 0.08
1.5	1000	64 ± 7	2.01 ± 0.06	1.55 ± 0.07	1.20 ± 0.09
1.5	2000	87 ± 12	2.06 ± 0.04	1.60 ± 0.07	1.22 ± 0.08
1.5	3000	117 ± 14	2.11 ± 0.04	1.65 ± 0.06	1.25 ± 0.08
1.5	4000	135 ± 21	2.15 ± 0.04	1.67 ± 0.05	1.26 ± 0.06
1.5	5000	159 ± 27	2.18 ± 0.03	1.69 ± 0.05	1.27 ± 0.07

A universal power law in the spatial distributions of pancreatic islets has been reported in guinea pigs, dogs, pigs, goats and cows with the fractal dimension D_2S = 1.5, measured in 2D pancreatic sections. In this study, we demonstrated that humans also followed the universal power law with the same fractal dimension. Furthermore, with the advanced imaging method, OPT, we found the fractal dimension in the 3D pancreas to be D₃ = 2.1 in mice. Note that the null hypothesis, D₃ = 2 (surface in the Euclidean space), was rejected (P < 0.001).

The pancreas has distinct developmental origins of the dorsal (splenic) and ventral (duodenal) pancreas, which later fuse [22]. The gastric lobe is formed by a perpendicular growth from the dorsal pancreas a few days after its formation [23]. Our recent study has shown regional differences in the pancreas that the gastric lobe has a higher relative number of islets (i.e. islets mm⁻³ of pancreatic tissue) than the splenic and duodenal lobes [24]. However, we found indistinguishable fractal dimensions of spatial islet distributions in all three lobes of the pancreas, although the higher islet density in the gastric lobe was reflected as a slightly higher fractal dimension that was not statistically significant. This suggests that the three regions may be governed by the same developmental rule.

The fractal dimension in the 3D pancreas has been extrapolated as D₃ = 2.5 by adding a unit co-dimension to the 2D section fractal dimension D_2S = 1.5, under the assumption of symmetric islet distributions in space [8]. Because the diffusion-limited aggregation (DLA) in 3D generates a tree-like structure with the fractal dimension D₃ = 2.4 [25], it has been proposed as the underlying mechanism of the spatial distribution of pancreatic islets [18]. However, our direct measurement of the 3D fractal dimension, D₃ = 2.1, was clearly different from D₃ = 2.4 of the DLA. Furthermore, endocrine cells are continuously nucleated from the source of the endodermal epithelium during pancreatic development. Therefore, the organ expansion is fundamentally different from the DLA in which particles scattered in a large volume aggregate to the center of a cluster [26]. The self-avoiding growth model, however, captured the organ expansion process. The DLA-like branching pattern has also been described by a neurite growth model attracted (not avoided) by chemical gradients [27]. Our model also explained the symmetry breaking, D₃ ≠ D_2S + 1, because the self-avoiding growth prefers to grow in directions as orthogonal as possible to the locations of the proximal existing nodes. When the self-avoiding effect becomes negligible (smaller m), spatial symmetry is restored.

The model with m = 1.5 could reproduce similar 3D, 2D projections and 2D section fractal dimensions of the spatial islet distributions in mice, considering the uncertainty of measurement. In particular, the uncertainty in the 2D section is considerable because the fractal dimension of the 2D section depends greatly on the position and orientation of the section plane in the pancreas volume. In general, smaller animals had smaller fractal dimensions (D_2S) in the spatial islet distributions: mice (1.31 ± 0.07), rats (1.42 ± 0.04) [9] and humans (1.53 ± 0.07). To reproduce the spatial islet distributions in larger animals including humans, the same model with a smaller m (weaker self-avoiding effect) could recreate a larger fractal dimension in the larger pancreas with more islets. This growth algorithm can be generally applicable to generate self-similar structures such as the vasculature with flexible fractal dimensions.

The inter-species differences of lung airway structures have been carefully examined by considering a multiplicity of scales [28]. Outer scales and, in more generality, the question of multi-fractal scaling would indeed be an interesting explanation for the islet data in different species. However, our theoretical attempts were limited by the amount of data we had available. The species difference can be simply understood in terms of self-avoidance, if we accept that different species can have different strengths of self-avoidance in islet development. Our data are limited to just center coordinates of islets, unlike the lung data with specific branch lengths and diameters. Nevertheless, multi-fractal scaling might lead to a model without self-avoidance that might apply to all species.

Node generation and self-avoidance in the phenomenological model could reflect the dual roles of blood vessels for inducing and inhibiting endocrine cell differentiation. Nevertheless, this model gives little insight into the specific mechanisms of self-avoidance. In addition, the model did not generate a conserved fractal dimension during its evolution. Ignoring the measurement uncertainty, the fractal dimension does not change with age in guinea pigs [9] and also in the mice shown in this study. During the evolution of the model, however, the fractal dimension slightly increased until the node number became large enough to overcome finite size effects.

The self-avoidance of islet formation may play a role for producing islets with optimal sizes. Endocrine cells tend to aggregate in embryonic development [29] and even in vitro culture [30]. The continuous generation of proximal cells, therefore, could lead to a gigantic cluster of cells. Constraints on endocrine cell differentiation by blood vessels [21] may inhibit the formation of large aggregates of cells and help to generate optimal sizes of islets evenly distributed throughout the pancreas. The functional importance of islet sizes has been emphasized [31–35]. As progress is made in both tissue scaffold engineering and the induction of desired lineage characteristics in induced pluripotent stem cells, the geometry of how functional replacement organs might be best engineered becomes important. As such, our growth model provides a simple characterization of how to grow finite samples with the desired geometric properties for the optimal function of pancreatic islets.

The work was supported by funding from the Kempe Foundations, Umeå University, and the Intramural Research Program of the National Institutes of Health, National Institute of Diabetes and Digestive and Kidney Diseases.