Reconstructing human pancreatic islet architectures using computational optimization

We outline a general methodology based on computational optimization and experimental data to reconstruct human pancreatic islet architectures. By using the nuclei coordinates of islet cells obtained through DAPI staining, cell types identified by immunostaining, and cell size distributions estimated from capacitance measurements, reconstructed islets composed of non-overlapping spherical cells were obtained through an iterative optimization procedure. In all cases, the reconstructed architectures included >99% of the experimental identified cells, each of them having a radius within the experimentally reported ranges. Given the wide use of mathematical modeling for the study of pancreatic cells, and recently, of cell-cell interactions within the pancreatic islets, the methodology here proposed, also capable of identifying cell-to-cell contacts, is aimed to provide with a framework for modeling and analyzing experimentally-based pancreatic islet architectures.


Introduction
Pancreatic islets constitute the endocrine part of the pancreas and are essential for glucose homeostasis. It has been estimated that a healthy human pancreas has ∼3 million islets dispersed throughout the pancreas with a mean diameter of 108.9 ± 6.2 μm. 1 Pancreatic islets are mainly composed of insulinproducing β-cells (∼65%), glucagon-producing αcells (∼30%) and somatostatin-producing δ-cells (∼5%), 2,3 with other cells such as the ghrelingproducing ε-cells 4 and the pancreatic polypeptideproducing (PP-) cells 5 also present but in a much lower proportion than the α, β and δ-cells. 6 Insulin is the only hormone capable of lowering glucose levels directly and glucagon is the main hyperglycemic hormone. 7 Somatostatin, on the other hand, inhibits both insulin and glucagon secretion, thus participating indirectly in the regulation of glucose levels. [8][9][10] At the single-cell level, pancreatic α, β and δ-cells rely on the generation of an electrical activity pattern for the secretion of glucagon, insulin and somatostatin, respectively, although it has been shown recently that Ca 2+ plays a relevant role in δ-cells independently of the electrical activity. 11 The stimulus-secreting processes in the α, β and δ-cells can be found in recent works on the subject (see for instance refs. 8,[12][13][14][15][16][17]. At a higher level of organization (i.e. islets), pancreatic cells communicate via direct electrical coupling through gap-junctions (between β cells and β and δ cells), 18,19 as well as paracrine, autocrine and juxtacrine signaling. [20][21][22][23] These interactions between cells depend on the number, distribution and organization of the α, β and δ-cells within the islets. 2,24,25 Several factors are involved in the formation, development and organization of pancreatic islets. For instance, it has been shown that factors such as the hepatocyte nuclear factor 6 26 and Nkx2.2, 27,28 as well as the mTOR pathway 29 and the Roundabout receptors 30 have an important role for the spatial organization of cells within the islet. It is also known that the architecture of pancreatic islets is disrupted in type 2 diabetes in rodents, non-human primates and humans, [31][32][33] including a decrease in β-cell mass in people living with type 2 diabetes 2,32,33 and the formation of amyloid deposits within the islets. [34][35][36] In addition to the functional implications of a loss of βcell mass (i.e. impaired insulin secretion), it is reasonable to hypothesize that alterations of the islet architecture also affect the interactions between cells within the pancreatic islets.
For more than three decades, mathematical modeling has been widely used to study both isolated pancreatic cells, 18,[37][38][39][40][41][42][43] clusters of β-cells [44][45][46][47][48][49][50] and more recently, models of pancreatic islets. 18,[51][52][53] While most recent models have been developed using real three-dimensional structures of both mice and human islets in terms of the nucleus position and identity of the cells (see for instance refs. 18,51,54 ) an aspect not explicitly considered in the models of pancreatic islets is the heterogeneity of sizes of the different cells. For this reason, in this work we develop a methodology based on computational optimization to generate human islet architectures using the distribution of radius of the α, β and δ-cells obtained from electrophysiological capacitance measurements, the position of the cells' nucleus obtained through DAPI staining, and the type of cells determined by immunostaining.

General methodology
The methodology proposed aims to reconstruct the architecture of pancreatic islets based entirely on experimental data. In summary, as depicted in Figure 1, an initial islet (iteration 0) is proposed by assigning a radius (r j,0 ) to each cell (j) following the radii distribution of the corresponding cell type, assuming that initially the cell nuclei, with coordinates (x n,j,0 , y n,j,0 , z n,j,0 ), are located at the center of the cells with coordinates (x c,j,0 , y c,j,0 , z c,j,0 ), where the index j identify the cells and i = 0 represents the iteration number. Then, the Simulated Annealing algorithm is used to find the islet with the minimum possible number of overlapped cells. As described graphically in Figure 3, at each iteration of the optimization algorithm a cell is randomly selected, a new radius is assigned, and its center is slightly displaced in the vicinity of the cell nucleus in a randomly selected direction (x, y or z) while keeping the nucleus coordinates fixed and always enclosed inside the cell (i.e. generating a new islet). This process is repeated until the predefined stop criteria are met and the optimized islet is then given as a result (details about the optimization process are given in Section 2.3). In the case of remaining pairs of overlapped cells at the end of the optimization procedure, an additional postprocessing step was performed (see Figure 4), consisting in randomly removing one cell from each pair of overlapped cells in order to obtain the final islet architecture.

Experimental data
In previous works, Hoang et al. 51,55 provided with six three-dimensional structures of human islets (fully available without restrictions under the terms of the Creative Commons Attribution License), consisting in the coordinates of each DAPI-stained nuclei along with the type of cell Figure 1. General methodology. An initial pancreatic islet is proposed based on nuclei coordinates obtained through DAPI staining (x n,j,i , y n,j,i , z n,j,i ), cells types determined by hormone staining (α,β,δ), and radii distributions obtained from capacitance measurements. Once the number of overlapped cells has been calculated for the initial islet N o (0), it is then used as the starting point for the Simulated Annealing algorithm in which, at each iteration (i), a random cell is selected and its radius and center coordinates are modified in order to generate a new islet. The number of overlapped cells at iteration i (N o (i)) is then calculated and the new islet is either accepted or rejected in accordance with the established criteria (see main text for details). Once a stop criteria is satisfied, an optimized islet is obtained. The final islet is obtained after a postprocessing step in which one of each pair of the remaining overlapped cells are randomly selected and neglected. obtained by immunostaining (a detailed description of the methodology used can be found in the original work by Hoang et al. 55 ) The main characteristics of the structures of human islets used in this work are summarized in Table 1.
In a recent work, Camunas-Soler et al. 56 measured the capacitance of 1369 human islet cells. By assuming sphericity, and using the widely known relation between membrane capacitance and surface area of 1 μF/cm 2 , the capacitance measurements by Camunas-Soler et al. 56 allowed us to estimate the size of the human α, β and δ-cells. Median values of 3.2, 5.3 and 3.9 pF were reported for the capacitances of the α, β and δ-cells, corresponding to spherical cells of radius 5, 6.5 and 5.6 μm, respectively. The corresponding minimum sizes reported for the human α, β and δ-cells were 4.3, 4.5 and 4.4 μm (2.4, 2.5 and 2.4 pF) while the maximum sizes were 6.9, 8.2 and 6.9 μm (5.9, 8.5 and 6 pF), respectively.

Optimization of islet architectures
The Simulated Annealing (SA) algorithm 57 was adapted to determine the cells' radii while minimizing the number of overlapped cells in the reconstructed islets. In short, the SA optimization algorithm is a process in which at each iteration a new set of parameters is randomly generated, accepting all the points that yield a lower value for the defined objective function, but also, accepting points that increase the objective function with certain probability depending on a parameter T (called "temperature" in the work by Corana et al. 57 ) as a way to prevent the process from getting trapped in a local minima, thus being more likely to reach the global minima.
As described schematically in Figure 2, the SA algorithm was implemented in such a way that the objective function to minimize is the number of overlapped cells (N o ), and the points randomly generated at each iteration i are the cells' radii (r j,i ) and the center coordinates (x c,j,i , y c,j,i , z c,j,i ).
An initial islet was generated using the experimental information available (i.e. the nucleus coordinates, cell types and cell size distributions). First, cells were centered at the nucleus coordinates and initial radii r j,0 were assigned randomly to each cell j in accordance with the corresponding radii distribution for the different cell types (α, β or δ). Radii were assigned by generating normal distributed random numbers with mean and standard deviation in accordance with the experimental reports, accepting only values laying within the range of values given by Camunas-Soler et al. 56 (see Section 2.2).  51 For the sake of comparison, the final characteristics of the reconstructed islets are also shown.

Experimental human islets
Final reconstructed islets δ  1  588  149  318  121  583  148  316  119  2  2264  430  1468  366  2252  427  1461  364  3  3253  1092  1543  618  3223  1082  1523  617  4  3544  969  2225  350  3516  961  2209  346  5  2096  649  1173  274  2084  642  1168  274  6  2858  837  1361  660  2841  830  1355  656 N: total number of cells; N ɑ , N β, N δ : number of ɑ, β and δ-cells, respectively At the initialization step, as well as at every iteration of the optimization process, the number of overlapped cells (i.e. the objective function, N o (i)) was calculated as: where i indicates the number of iteration (0 being the initial islet), N is the number of cells in the islet, r j,i and r k,i are the radius of j and k cells, and d j,k is the Euclidean distance between the corresponding center coordinates at iteration i. In practice, Eq. 1 means that when the sum of radii of the j and k cells is greater than the distance between their center coordinates, the number of overlapped cells, N o (i) increases by one. Thus, at each iteration i, Eq. 1 yields the total number of overlapped cells. As explained schematically in Figures 2 and 3, at each iteration of the SA algorithm a new islet was generated by randomly selecting a cell (j), assigning a new radius, r j,i+1 , in the range of the experimental reports for the type of cell selected, randomly selecting a direction of movement of the center coordinates q c,j,i (q = x, y, z), and calculating the displacement of the cell center randomly as: where U represents a uniform distribution in the given range. Moving the center of the selected cell according to Eq. 2 ensures that the nucleus remains inside the cell, while allowing the center of the cell to move in the vicinity of nucleus coordinates.
The new islet was evaluated by calculating the number of overlapped cells using Eq. 1. If that is, if the number of overlapped cells at the i + 1 iteration is lower than the number of overlapped cells at the previous iteration, the new islet is accepted. Otherwise, the new islet . Generation of new islets during the iterative optimization process. At each iteration i of the Simulated Annealing algorithm, a random cell j with radius r j,i (red dashed line), nucleus coordinates (x n,j,i , y n,j,i , z n,j,i ) shown in orange, and center coordinates (x c,j,i , y c,j,i , z c,j,i ) shown in yellow, is selected. First, a new radius r j,i+1 is determined, then the direction of movement of the center coordinates is selected and finally new center coordinates are generated. In the example shown, the center of the cell is displaced in the y direction.
Note that in Eq. 3 the acceptance probability decreases either as T decreases or as the difference A decrease of the parameter T is determined by a scheme consisting in proposing an initial value (T = 50) and then decreasing its value by a factor of 2 either when a predefined maximum number of islets (max acc ) was accepted or when a predefined maximum number of iterations (max iter ) was reached for each value of T. Values for max iter and max acc were defined in terms of the number N of cells in each islet (max acc = 500 × N and max iter = 1000 × N). This process (i.e. decreasing T and either accepting max acc islets or reaching max iter iterations) was performed iteratively until a predefined stop criteria was satisfied. Here, we established a stop criteria based on the maximum and minimum number of overlapped cells at iteration i (max{N o (i)} and min{N o (i)}, respectively) as: where � is a tolerance parameter with a value of 0.005. If the stop criteria is satisfied and N o �0, a postprocessing step was performed as described schematically in Figure 4. The postprocessing step simply consists in identifying the remaining overlapped pairs of cells and randomly selecting and removing one cell from each pair of overlapped cells. In the example schematically shown in Figure 4, three pairs of cells were overlapped at the end of the optimization procedure; therefore, in the postprocessing step one cell of each pair of overlapped cells were randomly selected and removed from the optimized islet.

Quantification of cell-to-cell contacts
Cell-to-cell contacts between α cells (N αα ), α and β cells (N αβ ), α and δ cells (N αδ ), β cells (N ββ ), β and δ cells (N βδ ) and δ cells (N δδ ) were quantified in the final optimized human islets as: with x = y = α, β, δ indicating the type of cell-to-cell contact to be quantified, r j,x j and r k,y k representing the radius of the j and k cells of type x j and y k respectively (α, β or δ), and s being a distance tolerance parameter determining the acceptable distance between cell surfaces to be considered as in contact (in this work we used a value of s = 1 μm). For instance, in order to quantify N αβ , the number of contacts between α and β-cells(i.e. x = α and y = β or y = α and x = β), N αβ would increase by one when the j and k cells are either α and β cells or β and α cells, respectively, and the sum of their radii plus the tolerance parameter s is equal or greater than the distance between the two cells.

Computational aspects
Code of the optimization process was written in C using the OpenMP library and run in the Yoltla cluster of the Laboratorio Nacional de

Results
Final reconstructed islets are presented in Figure 5, where the cell type, final radius, nucleus coordinates (black points), and center coordinates (red points) can be observed. Cross-sections of the six human islets reconstructed are shown in Figure 6 where it can be seen that cells do not overlap, which is one of the main advantages of the proposed methodology.
A summary of the optimization processes for the six human islets reconstructed are shown in Table 2. One important aspect shared by the six human islets obtained is that in all cases the islets included > 99% of the cell identified experimentally. The number of overlapped cells N o at the end of the optimization process ranged from 5 to 30 (see Table 2), being the mean distance between overlapped cells at the end of the optimization process between 2.49 and 3.65 µm, although the distance between nuclei of overlapped cells was as small as 0.18 µm.
The distributions of cells' radii for the six final islets are shown in Figure 7A and presented in detail in Table 3. Practically all the islets showed the same radius distribution for the three types of cells: α-cells with mean radii between 4.8 µm (islets 3-6) and 5.0 µm (islet 2), with minimum and maximum radii of 4.3 µm (all islets) and 6.8 µm (islets 1 and 3-6); β-cells with mean radii ranging between 5.1 µm (islets 5 and 6) and 5.4 µm (islet 1), with minimum and maximum radii of 4.5 µm (all islets) and 8.2 µm (islets 2-4 and 6), and δ-cells with mean radii in the range of 5 µm (islets 3, 4 and 6) and 5.2 µm (islet 1), with minimum and maximum radii of 4.4 µm (all islets) and 6.9 µm (islets 2-6). These distributions are in accordance with the experimental reports by Camunas-Soler et al. 56 (described in Section 2.2) who reported cells radii between 4.3 and 6.9 µm (α-cells), 4.5 and 8.2 µm (β-cells) and 5.9 and 6.9 µm (δ-cells). Also in Table 3, the total islet volumes, as well as the volumes of α, β and δ cells are shown. Ionescu et al. 1 found a mean islet volume of 686,994 ± 107,297 µm 3 , while the total volume of the reconstructed islets ranged between 374,747 and 2,028,662 µm 3 . It should be noted, however, that the human islets analyzed by Ionescu et al. had a mean diameter of 113.8 ± 7 µm, while Hoang et al. 51 described larger islets, as can be seen in Figure 5.
As expected, the computing time of the optimization procedure highly depended on the number of cells composing the islets. This can be seen when comparing the computing times of the smallest islet (islet 1) initially composed of 588 cells, and that of islet 4, initially composed of 3544 cells. While the computing time of the former was approximately 2 hours and 15 minutes, the optimization of the latter lasted more than 5 days and 19 hours. In all cases, several million iterations were performed (11.6 million for islet 1 and more than 50 million for islets 2-6, see Table 2).  Table 2). α, β and δ-cells are shown in red, green and blue respectively. The purple and red dots indicate the location of the center and nucleus coordinates, respectively.
Convergence plots of the optimization processes for the six islets reconstructed are shown in Figure 7B, where it can be seen how the number of overlapped cells decreased as the parameter T was reduced, which in turn reduced the probability of accepting islets that increased the number of overlapped cells.  Cell-to-cell contacts were determined for the six human islets reconstructed. As described in Table 4, total contacts between all types of cells ranged from 750 for islet 1, composed of 583 cells, to 5533 for islet 4, composed of 3516 cells. Homotypic contacts were the most abundant within the islets (57.1%-70.4%), with heterotypic contacts ranging between 29.6% and 42.9%. Homotypic contacts were predominantly between β-cells while the heterotypic contacts were mostly between α and β-cells. The number and percentages of all the cell-to-cell contacts can be consulted in Table 4. A visual description of these results can be seen in Figure 8, where contacts identified in islet 1 are shown (similar results obtained for the other five islets can be consulted in the repository related to this article). These results are in accordance with the reports by Hoang et al., 55 who found that the relative attraction between homotypic cells is significantly stronger than the attraction between heterotypic cells.

Discussion
In this work, we proposed a methodology based on computational optimization to reconstruct human islet architectures from experimental data. Even though we have focused on human islet architectures, the proposed methodology can be readily used to reconstruct islet architectures from other species. The resulting islet architectures, composed of nonoverlapping cells, reproduced key experimental observations such as the radii distribution and the location of the cells' nuclei. In addition, the reconstructed architecture was also used to obtain information about the cell-to-cell contacts, a key aspect for the study of the interactions between cells within the islets.  Camunas-Soler et al. 56 reported cell capacitances equivalent to minimum and maximum cell sizes o 4.3 and 6.9 μm for the ɑ-cells, 4.5 and 8.2 μm for the β-cells and 4.4 and 6.9 μm for the δ-cells. V tot is the total volume, V ɑ , V β , and V δ are the ɑ, β and δ-cells volumes. All volumes are given in μm 3 . Radii are given in μm. Table 4. Cell-to-cell contacts in the six human islets reconstructed. Throughout the years, several multicellular computational models of pancreatic cells have been proposed, aiming to capture important functional aspects related to the inter-islet communication between cells in the pancreatic islets. Even early models such as those developed by Chay and Kang, 44 Sherman, Rinzel and Keizer 45 and Sherman and Rinzel 46 have focused on capturing the effects of electrical communication between neighbor β-cells using regular arrays of same-size spherical cells. However, these early models, and even more recent models have not captured accurately the heterogeneity of cell sizes inherent of intact pancreatic islets. In spite of this, computational models of human and mice islets have been developed to give a better understanding of how electrical activity in islet cells is affected by gapjunction heterogeneities in health and disease, 47-49,58-60 the role of β-cell hubs in islet function, 54,61,62 the effect of metabolic defects in subpopulations of β-cells, 50 hormone pulsatility and the effects of cells organization in normal and diabetic islets 51 and how the coupling between βcells and β and δ-cells, regulate the activity of αcells via paracrine mechanisms. 18 Most of these models adopted simple geometrical representations of the clusters of cells such as cubes of same-size cells. Only recent models by Hoang et al., 51 Briant et al. 18 and Lei et al. 54 used real three-dimensional structures of mouse and human islets.
On the other hand, other models of β-cells and pancreatic islets based on the finite element method (FEM) have focused on simulating the spatial distribution of intra and extracellular signals. 39,52,53,[63][64][65] For instance, Buchwald et al. 52,53,63,64 developed spatial models of in vitro spherical islets to simulate oxygen and glucose consumption and insulin secretion using a simplified geometrical representation of the pancreatic islets, not including an explicit geometrical representation of the islet cells. 52,53,63 Although all these works based on simplified islet geometries yield very interesting results, we believe that by adopting a more realistic representation of islet cells will improve the reach of future theoretical work in the field of islet physiology.
In spite of the satisfactory results described throughout the article, there are some limitations that must be highlighted. Firstly, the space occupied by the vasculature and other structures has not been considered. Secondly, we have used capacitance measurements of dispersed cells to estimate the radii of the α, β and δ-cells, which could not be reflecting their real size in intact islets. Finally, it has been assumed that α, β and δcells can be represented by spheres. Although this might be a good approximation for α and β-cells,  Figures 5 and 6. Note that in order to highlight he cell-to-cell contacts, the cell's radii are not shown. α, β and δ-cells are shown in red, green and blue, respectively. Similar results for islets 2-5 can be consulted in the repository associated to this article (see Section 2.5).
which are rounded or rhomboid in shape, 8 it is known that some δ-cells show long neurite-like processes 8,66 that cannot be directly generated by the methodology here proposed. Despite this, the reconstructed islets had a total volume in the range of the experimental reports, 1 and given that the estimated cells' radii are based on capacitance measurements, the spherical δ-cells in the reconstructed islets (and also α and β-cells) had surfaces areas in accordance with the experimental observations.
Overall, given the increasing importance of modeling for the study of inter-islet communication within the islets of Langerhans, we believe that the methodology here proposed will be a useful tool for the development of more realistic models of pancreatic islets, suitable not only for models focused on the electrical behavior of coupled cells, but also for models of islets aiming to describe the spatial distribution of ions, hormones or other diffusible molecules. Furthermore, in the light of recent studies linking functional and molecular phenotypes at the single-cell level, 56 we believe that future computational models of pancreatic islets will be able to capture both the morphological and functional heterogeneities of islet cells in health and disease.