Computational Reconstruction of Pancreatic Islets as a Tool for Structural and Functional Analysis

Summary

In this protocol, the pancreatic islets are reconstructed and analyzed using computational algorithms implemented in a dedicated multiplatform application.

Abstract

Structural properties of pancreatic islets are key for the functional response of insulin, glucagon, and somatostatin-secreting cells, due to their implications in intraislet communication via electric, paracrine, and autocrine signaling. In this protocol, the three-dimensional architecture of a pancreatic islet is firstly reconstructed from experimental data using a novel computational algorithm. Next, the morphological and connectivity properties of the reconstructed islet, such as the number and percentages of the different type of cells, cellular volume, and cell-to-cell contacts, are obtained. Then, network theory is used to describe the connectivity properties of the islet through network-derived metrics such as average degree, clustering coefficient, density, diameter, and efficiency. Finally, all these properties are functionally evaluated through computational simulations using a model of coupled oscillators. Overall, here we describe a step-by-step workflow, implemented in IsletLab, a multiplatform application developed specifically for the study and simulation of pancreatic islets, to apply a novel computational methodology to characterize and analyze pancreatic islets as a complement to the experimental work.

Introduction

The pancreas is divided into regions referred to as head, neck, body, and tail, each having different structures, functions, and anatomical position^1,2. From a functional viewpoint, the pancreas can be divided into endocrine and exocrine systems with the former responsible for the secretion of hormones critically involved in the regulation of glucose homeostasis, while the latter contributes to food digestion via the secretion of enzymes into the duodenum¹. Pancreatic islets constitute the endocrine tissue of the pancreas and are responsible for the secretion of glucagon, insulin, and somatostatin, secreted from ɑ, β, and δ-cells, respectively³. In addition to their intrinsic regulatory mechanisms, these cells are regulated via direct electrical communication (between β-cells and likely β and δ-cells), and also by paracrine and autocrine signaling^4,5,6. Both mechanisms are highly dependent on the islet architecture (i.e., the composition and organization of the different types of cells within the islet)^7,8. Importantly, islet architecture is altered in the presence of diabetes, most likely disturbing intraislet communication as a result^9,10.

The study of pancreatic islets involves a wide range of experimental methodologies. Among these, the use of fluorescence techniques to determine the number, location, and type of the different cells in the islet has allowed to study the structural and morphological properties of pancreatic islets^11,12,13 and to gain a better understanding of the functional implications in health and disease. As a complement, computational models of pancreatic cells^14,15,16 and, more recently, pancreatic islets^12,17,18,19 have been used in the last decades to evaluate aspects difficult or even impossible to address experimentally.

In this protocol, we aim to bridge the gap between the experimental and computational work by outlining a methodology to reconstruct islet architectures, to analyze their morphological and connectivity properties through quantitative metrics, and to perform basic simulations to evaluate the functional implications of the islet properties.

The protocol described below is based on computational algorithms specifically designed for the study of pancreatic islets. In summary, in the first step of the protocol, the islet architecture is reconstructed from experimental data using the algorithm recently proposed by Félix-Martínez et al.¹⁹ in which nuclear positions obtained through 4′,6-diamidino-2-phenylindole (DAPI) staining and cellular types identified through immunofluorescence (as described in detail by Hoang et al.^11,12) are processed in an iterative optimization procedure. This leads to determining the optimal size and position of each cell and obtaining an islet composed of non-overlapping cells. Secondly, based on the reconstructed architecture, cell-to-cell contacts are identified to determine the connectivity properties and to generate the corresponding islet network which allows the user to obtain quantitative metrics to further describe the islet architecture (details about the reconstruction algorithm can be consulted in the original work on the subject¹⁹). Finally, basic functional simulations are performed using the modeling approach proposed by Hoang et al.¹² in which, based on the pulsatile nature of hormone secretion observed experimentally^20,21, each cell is treated as an oscillator, and therefore the islet is represented as a network of coupled oscillators following the connectivity properties of the reconstructed islet.

Given the computational complexity of the algorithms used in this protocol, all the steps involved have been implemented in a standalone application²² with the main objective of approaching these computational tools to all the interested readers regardless of their level of experience in the use of specialized software or programming languages.

Protocol

NOTE: A schematic diagram of the protocol is shown in Figure 1. A step-by-step description is given is as follows (see Supplementary File 1 for details about the control panels used at every step of the protocol).

Figure 1: Flow diagram. A flow diagram describing the sequential order of the protocol as implemented in IsletLab. Please click here to view a larger version of this figure.

1. Installing IsletLab in Linux

NOTE: Follow the instructions given in sections 2 and 3 of the Supplementary File 2 to install IsletLab in windows or macOS.

1. Open an internet browser and go to https://github.com/gjfelix/IsletLab. Download and extract the IsletLab repository files.
   NOTE: The Supplementary File 3 includes a copy of the version of IsletLab used in this protocol.
2. Verify that the gcc and nvcc compilers are installed. Open a terminal and type the following commands:
   gcc --version
   nvcc --version
   Follow the instructions listed in section 1 of the Supplementary File 2 if any of these commands is not recognized by the system.
3. Download and install the data science platform (see Table of Materials). Open a terminal and go to the IsletLab folder.
4. Create a new environment by typing the following command in the terminal:
   conda env create -f isletlabgui_v1.0.yml
5. Activate the new environment by typing:
   conda activate isletlab_v1.0
6. Launch the IsletLab application by typing the following command in the terminal (see Figure 2 for a description of the main window):
   python isletlabgui_v1.0.py

Figure 2: The user interface of IsletLab. The interface is composed of three main panels: configuration (1), statistics (2), and graphics (3) panels. The graphics toolbar (4) is located at the bottom of the graphics panel. Please click here to view a larger version of this figure.

2. Islet reconstruction

1. Prepare the input data. Organize the input islet data (i.e., cell types and three-dimension (3D) coordinates) in a four-column file in which column 1 contains the cell type (coded as 11: ɑ-cell, 12: β-cell, 13: δ-cell) and columns 2 to 4 contain the x, y, and z coordinates, respectively (see an excerpt of an input file in section 5 of the Supplementary File 2 or the input test file included in the IsletLab repository).
2. Ensure that the input file does not include column headers. Use the example input file included in the IsletLab repository to follow the protocol if needed (Input_Islet_test.txt).
   NOTE: Input data used to develop the algorithms implemented in IsletLab was obtained experimentally as described in detail by Hoang et al.^11,12 using DAPI staining to determine the nuclear positions and immunofluorescence to identify the cellular types.
3. Click the Load Initial Islet button and select the file containing input data to generate an initial islet, the 3D representation, and the corresponding statistics. In case of error, the following message will be displayed: 'Error loading islet file'. To correct this, repeat step 2.1.
4. Configure the reconstruction process. Click the Reconstruction Settings button and modify the optimization parameters (see sections 4 and 5 of Supplementary File 2 for a brief description of the algorithm and the parameters shown in the reconstruction settings window).
5. Set the initial temperature = 1, iterations factor = 1, and acceptance factor = 1. Click the OK button to save the parameter values.
   NOTE: As a rule of thumb, an increase in the temperature parameter, as well as in the iterations and acceptance factors, will yield better reconstruction results in terms of the experimental cells included in the reconstructed islet at the expense of an increase in the computing time.
6. Click the Reconstruct Islet button to open the reconstruction log window (see section 6 of Supplementary File 2 for a description of the information provided in the reconstruction log).
7. Click the Run button to start the reconstruction process. Monitor the optimization process until the message: 'Please close this window to continue' is displayed. Close the reconstruction log window to generate the visual representation of the reconstructed islet and to calculate the related statistics (see the left column in Figure 3A-C).
8. Evaluate the results of the reconstruction process by analyzing the optimization statistics shown in the final islet tab of the statistics panel. Particularly, focus on maximizing the percentage of experimental cells included in the reconstructed islets (% of experimental) or equivalently, on minimizing the number of overlaps (also shown graphically in the convergence plot, right column in Figure 3A-C).
9. If the % of experimental statistic is considered low according to the user objectives, restart IsletLab as described in step 9.1, and increase the initial temperature, iterations factor and acceptance factor in the reconstruction settings and repeat steps 2.1-2.4 until satisfactory results are obtained (> 95% of experimental cells). See sections 5-7 of the Supplementary File 2 for further description of the reconstruction process and the related results.

3. Identifying cell-to-cell contacts

1. Click the Reconstruction Settings button and set the Contact Tolerance parameter to define the cell-to-cell contact tolerance and click OK to save the parameter values.
   NOTE: The contact tolerance parameter represents the maximal distance between cells in order to be considered as in contact.
2. Click the Cell-to-Cell Contacts button to identify the cells in close contact. In the contacts tab, check for cells in contact shown graphically (black lines) in the graphics panel and the corresponding statistics (i.e., total, homotypic and heterotypic contacts and ɑ-ɑ, β-β, δ-δ, ɑ-β, ɑ-δ, β-δ contacts) shown in the statistics panel (see the left column of Figure 4A-C).
3. Increase (decrease) the contact tolerance parameter to increase (decrease) the number of cell-to-cell contacts. See section 8 of Supplementary File 2 for details about the identification of cell-to-cell contacts and the associated results.

4. Building the islet network

1. Click the Build Network button to generate the islet network and to calculate the associated network metrics (see the right column in Figure 4A-C).
   ​NOTE: A further description of the results related to the islet network generated is given in section 9 of Supplementary File 2.

5. Functional simulation of the reconstructed islet

1. Switch to the Simulation Tab of the configuration panel of the interface (Figure 5).
2. Select the desired mode of intrinsic frequency - Constant or Random - and click the Configure Intrinsic Frequency button to define the oscillators frequency (in Hz). If a random intrinsic frequency is selected, define the mean and deviation (in Hz) to generate normally distributed random frequencies (see section 11 of the Supplementary File 2 for a description of the simulation parameters).
3. Select the desired mode of the initial phase - Constant or Random. If a Constant initial phase is selected, click the Configure Initial Phase button to define the phase of the oscillators (in radians). If random initial phase is selected, system will assign random phases between 0 and 2π to all the oscillators.
   NOTE: If a constant initial phase is selected, all the oscillators will be initialized in phase.
4. Click the Configure Interactions button to define the cell-to-cell interaction parameters in the interaction strength window. For a description of the interaction parameters, see section 11 of Supplementary File 2.
5. Configure the simulation by defining the total simulation time (in s), time step (in s), and save factor (numbers of steps between saved data points). In the simulation shown in Figure 5 the total time was 20000 s, with a time step of 0.1 s and a save factor of 500.
6. Define the number of blocks, threads, and computing platform capability available to perform the simulation. Define these parameters in accordance with the specific characteristics of the graphics processing unit (GPU) available.
   NOTE: In the example shown in Figure 5, 36 blocks and 64 threads per block were used (2304 computing platform cores) since a GPU with 36 multiprocessors and 64 computing platform cores per multiprocessors was used. The value of the computing platform capability parameter used was 75 since the hardware computing platform capability was 7.5 (see section 10 of Supplementary File 2 for details about these parameters).
7. Click the Run Simulation button to open the simulation log window. Click the Run button to start the simulation and monitor the process until the legend: 'Please close the window to continue' is displayed. See section 12 of Supplementary File 2 for details about the simulation log.
8. Close the simulation log window to observe the simulation results (see Figure 5). Details about the simulation results are given in section 13 of Supplementary File 2.

6. Save the project (optional)

1. Click File > Export Project in the menu bar. Select the directory in which the project file will be saved and click the OK button.
   NOTE: The name of the project file is determined automatically based on the name of the initial data file. If the project is not saved, all the results and related files will be deleted automatically.
2. Load an exported project by clicking File > Load Project.

7. Save figures (optional)

1. Click the Save Plot Icon located in the plot toolbar to save the current visualization to an image file. Repeat this step for all the figures created throughout the protocol.
   ​NOTE: Plots and islet visualizations can be modified using the icons available at the plot toolbar of the graphics panel.

8. Load project (optional)

1. Click File > Load Project in the main menu and select a project file saved previously as described in steps 6.1-6.2.
   ​NOTE: A project will not be loaded correctly if the project file has been modified externally. If a project is loaded successfully, only the simulation capabilities will be available.

9. Restart the reconstruction process and analysis (optional)

1. At any time the user can select File > Restart to discard all the results and start a new reconstruction and analysis.

Representative Results

 The reconstruction of pancreatic islets using the methodology proposed by Félix-Martínez et al.¹⁹ is highly dependent on the parameters given to the optimization algorithm (defined in the reconstruction settings). An example of this is shown visually in Figure 3 where reconstructed islets obtained using different sets of parameters are shown. First, in Figure 3A, a reconstruction that included 86.6% of the cells included in the initial data is shown (509 out of 588 cells, initial temperature = 1, iterations factor = 1, acceptance factor = 1). When the initial temperature and the iteration and acceptance factors are increased (initial temperature = 10, iterations factor = 5, and acceptance factor = 5, Figure 3B), a higher percentage (93.37%) of initial cells were included in the reconstructed islets (i.e., 549 out of 588 cells). Even better results can be obtained if much higher values are used, particularly for the iteration and acceptance factors, as illustrated in Figure 3C (initial temperature = 10, iterations factor = 1000, acceptance factor = 500), where the reconstructed islet is composed of 99.15% of the initial cells (583 out of 588 cells). The convergence plots (right column in Figure 3A-C), showing the evolution of overlapped cells as a function of the temperature, must be evaluated to determine how the parameters affect the optimization process. As a rule of thumb, the interaction and acceptance factors must be increased when the reconstructed islet includes a low percentage of initial cells. Consequently, the computing time will inevitably increase, since these factors directly increase the number of iterations evaluated. For instance, the computing time of the first reconstruction described above was 6 s. In contrast, the computing times of the second and third reconstructions were 21 s and 24 min 6 s, respectively.

Figure 3: Islet reconstruction using suboptimal sets of parameters in the reconstruction settings Using suboptimal sets of parameters could lead to a low percentage of experimental cells in the reconstructed islets. (A) Left: 86.6% of experimental cells were included in the reconstructed islet (initial temperature = 1, iterations factor = 1, acceptance factor = 1, computing time = 6 s). Right: convergence plot of the reconstruction process. (B) Left: 93.4% of experimental cells were included in the reconstructed islet (initial temperature = 10, iterations factor = 10, acceptance factor = 5, computing time = 21 s). Right: convergence plot of the reconstruction process. (C) Left: 99.15% of experimental cells were included in the reconstructed islet (initial temperature = 10, iterations factor = 1000, acceptance factor = 500, computing time = 24 min, 8 s). Right: convergence plot of the reconstruction process. The arrows in the convergence plots indicate the initial and final number overlapped cells of the reconstruction process (before the postprocessing phase of the reconstruction algorithm). Please click here to view a larger version of this figure.

Identification of cell-to-cell contacts from the reconstructed islet depends on the value of the contact tolerance parameter (defined in the reconstruction settings) as illustrated in Figure 4A-C, where the cell-to-cell contacts (represented by black lines), identified from the reconstructed architectures shown in Figure 3A-C, are presented. For instance, if a contact tolerance of 1 µm is defined, as in Figure 4A, only 290 cell-to-cell contacts are identified. In contrast, if the contact tolerance is increased to 2 µm, as in Figure 4B,C, the total contacts identified increased to 636 and 731, respectively (see the statistics panel in Figure 4A-C). These differences can also be noticed in the visual representation of the cell-to-cell contacts shown in the left column of Figure 4A-C, as the number of contacts between cells clearly increases as a higher value of the contact tolerance is used. It is worth highlighting that the number of contacts also depends on the number of cells included in the reconstructed islets, and therefore, the combination of the temperature parameter, iterations and acceptance factors, and contact tolerance ultimately determine the connectivity of the reconstructed islet, which is reflected on the islet networks formed and the corresponding network metrics, as shown in the right column of Figure 4A-C. The network plot allows the user to visualize how the different cells are connected. Quantitatively, connectivity properties of the islet are described in terms of the following network metrics: average degree, density, average clustering coefficient, efficiency, and diameter (details about these metrics can be consulted in section 9 of the Supplementary File 2).

Figure 4: Effect of the contact tolerance parameter in the identification of cell-to-cell contacts. (A-C) Left: cell-to-cell contacts identified from the reconstructed islets shown in Figure 3A-C (290, 636, and 731 total contacts in panels A, B, and C, respectively). Values used for the contact tolerance parameter were 1 µm (A) and 2 µm (B and C). Note that the number of cells included in the reconstructed islets also affects the number of cell-to-cell contacts identified. Right: networks generated from the cell-to-cell contacts are shown in the corresponding left column. Note the impact of the connectivity on the network metrics is highlighted in the statistics panel. Please click here to view a larger version of this figure.

Finally, once the islet has been reconstructed and the cell-to-cell contacts have been identified, a functional simulation can be performed (only when a compatible GPU is available). Typical simulation results are shown in Figure 5, including the summed oscillations of the different cell populations (ɑ, β, and δ-cells) and the whole islet (upper plot of the graphics panel in Figure 5). This figure shows the phase differences over time between the different cell populations as a result of the connectivity and interaction properties and allows the user to determine the contribution of each cell population (red, green, and blue lines) to the oscillatory behavior of the whole islet (black line). For instance, the upper panel of Figure 5 suggests that, at the population level, ɑ and β-cells oscillate completely out of phase, while δ-cells oscillate out of phase with ɑ and β-cells. Moreover, according to the simulation, the oscillatory behavior of the islet is dominated by the oscillations of the ɑ-cells, although the effect of the other cell populations can be also noticed. Note that the oscillatory signals of all the islet cells are saved automatically in a data file (see Table 1 and section 13 in Supplementary File 2), thus allowing the user to perform a detailed analysis of the simulation results. As a complement, the islet synchronization index, which reflects the phase coherence of the oscillations, is also calculated and displayed (bottom plot of the graphics panel in Figure 5). Note that the synchronization index ranges from 0 to 1, where 0 and 1 indicate a null and total synchronization between all the cells in the islet, respectively. The synchronization index plot can be therefore interpreted as a visualization of how the synchronicity between islet cells varies over time as a result of the connectivity and interaction properties of the reconstructed islet. Since the simulation performed is based on the idea of coupled oscillators¹² and heavily depends on the connectivity of the reconstructed islet, it is key to reach an acceptable islet reconstruction and cell-to-cell connectivity before performing a functional simulation.

Figure 5: The simulation parameters are defined in the configuration panel of the simulation tab. The results of the simulation are shown in the simulation tab of the graphics panel where the summed oscillatory behavior of the different populations of cells (ɑ, β, and δ) and the whole islet are shown (top). The synchronization index, a measure of phase coherence between the islet cells, is also shown (bottom). Please click here to view a larger version of this figure.

It is worth mentioning that practically at every step of the process, data files are generated. A description of the data files generated can be found in Table 1 and throughout Supplementary File 2.

File	Description
IsletFileName	Input data (given by the user)
IsletFileName_initial.txt	Initial islet architecture proposed by the algorithm to as initial step of the reconstruction
IsletFileName_reconstructed.txt	Reconstructed islet (not posptrocessed)
IsletFileName_postprocessed_islet.txt	Final reconstructed islet and postrocessed islet
IsletFileName_processlog.txt	Reconstruction log (optimization algorithm)
IsletFileName_overlapped_cells.txt	Overlapped cells at the end of the reconstruction process (postprocessing)
IsletFileName_all_contacts.txt	Adjancency matrix of all contacts
IsletFileName_aa_contacts.txt	Adjancency matrix of ɑ-ɑ contacts
IsletFileName_ab_contacts.txt	Adjancency matrix of ɑ-β contacts
IsletFileName_ad_contacts.txt	Adjancency matrix of ɑ-δ contacts
IsletFileName_bbbd_contacts.txt	Adjancency matrix of β-β and β-δ contacts
IsletFileName_bb_contacts.txt	Adjancency matrix of β-β contacts
IsletFileName_bd_contacts.txt	Adjancency matrix of β-δ contacts
IsletFileName_dd_contacts.txt	Adjancency matrix of δ-δ contacts
IsletFileName_Kmat.txt	Interaction matrix used in the simulation
IsletFileName_kuramoto_angles.txt	Results of the Kuramoto simulation

Table 1: Description of files saved as a part of the project file. Note that the file name used to save the project files are automatically defined by the initial data file selected by the user.

Supplementary File 1: Graphical description of the protocol using the control panels of IsletLab. Please click here to download this File.

Supplementary File 2: IsletLab documentation. Please click here to download this File.

Supplementary File 3: Includes all the files needed to install IsletLab. Please click here to download this File.

Discussion

The above protocol outlines a practical approach to reconstruct and analyze pancreatic islet architectures using novel computational algorithms. The main objective of this work is to enable the islet research community to derive quantitative metrics to characterize the morphological and connectivity properties of pancreatic islet architectures and to evaluate the possible functional implications of such properties via computational simulations.

While the algorithms adopted in this protocol have been previously described in detail^12,19, a direct and user-friendly implementation was lacking due to their relative complexity, thus limiting their use as a complementary tool to the experimental and theoretical work.

Firstly, a recent algorithm proposed by Félix-Martínez et al.¹⁹ is used to reconstruct the islet architectures from experimental data (e.g., nuclear coordinates and cell type). As a result, the user obtains an islet architecture composed of non-overlapping spherical cells with radii automatically assigned in accordance with the reported experimental distributions. In practice, the reconstruction algorithm is an iterative optimization procedure that becomes expensive from the computational standpoint as the number of cells in the islet increases. For this reason, it is highly recommended to use a multiprocessor system to take advantage of the parallel processing implementation of the algorithm described in this protocol. As described above, a key step for the reconstructing process is to define appropriate values for the parameters involved (i.e., iterations, acceptance factors, and initial temperature), since the computing time will be directly related to the number of iterations performed, in addition to the number of parallel processes used (i.e., threads parameter in the reconstruction settings). If computing time is not an issue, we strongly recommend using the highest values possible for the iterations and acceptance factors in order to increase the number of iterations performed.

The next steps of the protocol are the identification of cell-to-cell contacts and the generation of the islet network. Both steps are directly related to the reconstruction process and as such, the number of cells included in the reconstructed islets (and therefore the parameters involved), as well as the value of the contact tolerance used, is key to obtaining the best results possible.

Finally, if desired by the user, functional simulations can be performed through the implementation of the model of coupled oscillators proposed by Hoang et al.¹² using the connectivity network derived from the reconstruction process to configure the oscillatory system. Given that the simulation process involves solving a system of hundreds or thousands of coupled differential equations (one for each cell in the islet), the simulation algorithm has been implemented taking advantage of the possibility of performing parallel computations using the GPU, thus allowing the user to simulate considerable long simulations in a relatively short computing time. Key steps in the simulation stage of the protocol are to determine the appropriate number of blocks and threads available in the computing platform settings section of the simulation panel, an aspect directly related to the characteristics of the hardware used. The other parameters involved (intrinsic frequency, initial phase, and interaction strengths in the simulation panel), although relevant for the simulation results, are mainly related to the problem under investigation and must be defined by the user after thoughtful consideration in order to represent the desired simulation scenario.

Despite the advantages offered by the protocol, some limitations must be acknowledged. Firstly, the parameters related to the reconstruction process and identification of the cell-to-cell contacts are not unique and might vary from case to case. For this reason, although a rule of thumb can be used to determine the value of the required parameters, a trial-and-error approach is still unavoidable. Another aspect that could limit the applicability of the protocol is the computational resources needed, particularly for the reconstruction and simulation stages of the protocol. Despite these limitations, the fact that programming knowledge is not needed for the implementation of the protocol allows researchers from diverse backgrounds to readily make use of the proposed algorithms that otherwise would remain obscure for the non-specialized user.

Potential uses of the proposed protocol include the visualization of experimental data, comparative analysis of normal and altered islets (e.g., in the presence of type 1 or 2 diabetes), or even comparison between islets from different species using quantitative morphological, structural, and network-based metrics²³. Moreover, reconstructed islets using the protocol outlined here can be readily used to generate detailed functional mathematical models in which the connectivity and cell sizes determined by the reconstruction algorithm are complemented with detailed electrophysiological models of pancreatic cells to elucidate the functional implications of the intercellular communication within reconstructed islets.

Materials

Name	Company	Catalog Number	Comments
CUDA-capable NVIDIA graphics card			Required for the functional simulations
IsletLab			https://github.com/gjfelix/IsletLab (Follow the instructions to download and install the application.)