Experimental evaluation and computational modeling of the effects of encapsulation on the time-profile of glucose-stimulated insulin release of pancreatic islets

Background In type 1 diabetic patients, who have lost their ability to produce insulin, transplantation of pancreatic islet cells can normalize metabolic control in a manner that is not achievable with exogenous insulin. To be successful, this procedure has to address the problems caused by the immune and autoimmune responses to the graft. Islet encapsulation using various techniques and materials has been and is being extensively explored as a possible approach. Within this framework, it is of considerable interest to characterize the effect encapsulation has on the insulin response of pancreatic islets. Methods To improve our ability to quantitatively describe the glucose-stimulated insulin release (GSIR) of pancreatic islets in general and of micro-encapsulated islets in particular, we performed dynamic perifusion experiments with frequent sampling. We used unencapsulated and microencapsulated murine islets in parallel and fitted the results with a complex local concentration-based finite element method (FEM) computational model. Results The high-resolution dynamic perifusion experiments allowed good characterization of the first-phase and second-phase insulin secretion, and we observed a slightly delayed and blunted first-phase insulin response for microencapsulated islets when compared to free islets. Insulin secretion profiles of both free and encapsulated islets could be fitted well by a COMSOL Multiphysics model that couples hormone secretion and nutrient consumption kinetics with diffusive and convective transport. This model, which was further validated and calibrated here, can be used for arbitrary geometries and glucose stimulation sequences and is well suited for the quantitative characterization of the insulin response of cultured, perifused, transplanted, or encapsulated islets. Conclusions The present high-resolution GSIR experiments allowed for direct characterization of the effect microencapsulation has on the time-profile of insulin secretion. The multiphysics model, further validated here with the help of these experimental results, can be used to increase our understanding of the challenges that have to be faced in the design of bioartificial pancreas-type devices and to advance their further optimization. Electronic supplementary material The online version of this article (doi:10.1186/s12938-015-0021-9) contains supplementary material, which is available to authorized users.

is briefly summarized in the Method section; a fully detailed description is presented here.

Reaction rates
For all species, reaction rates (i.e., consumption and release rates, R) are assumed to follow Hilltype dependence on the local concentrations (i.e., generalized Michaelis-Menten kinetics): This functional dependence uses three parameters: R max , the maximum reaction rate [mol⋅m -3 ⋅s -1 ], C Hf , the concentration corresponding to half-maximal response [mol⋅m -3 ], and n, the Hill slope characterizing the shape of the response. This function (originally introduced by A. V. Hill [5,6]) is widely used for biological / pharmacological applications [7] as it provides transition from zero to a limited maximum rate via a smooth an continuously derivable function of adjustable width. The well-known two-parameter Michaelis-Menten equation [ perifusion studies: one concentrating on the effect of glucose using isolated human islets [9] and one concentrating on the effect of hypoxia using isolated rat islets [10].

Oxygen consumption
For oxygen consumption, the basic parameter values are maintained from our previous models [1,11]: n oxy = 1, R max,oxy = -0.034 mol⋅m -3 ⋅s -1 , and C Hf,oxy = 1 μM (corresponding to a partial oxygen pressure of p Hf,oxy = 0.7 mmHg) ( Table 1). By all indications, the assumption of a regular Michaelis-Menten kinetics (i.e., n oxy = 1) gives an adequate fit. Accordingly, at low oxygen concentrations, where cells only try to survive, oxygen consumption scales with the available concentration c oxy and, at sufficiently high concentration, it plateaus at a maximum (R max ). To account for the increased metabolic demand of insulin release and production at higher glucose concentrations, a dependence of R oxy on the local glucose concentration is assumed via a modulating function ϕ o,g (c gluc ): This is done to accommodate experimental indications showing increased oxygen consumption rate in islets when going from low to high glucose (see [1] and references therein). It is assumed that the oxygen consumption rate contains a base-rate and an additional component that increases due to the increasing metabolic demand in parallel with the insulin secretion rate (see eq. A6 later) as a function of the glucose concentration: In the present model, the base and metabolic components are assumed to be equal (ϕ base = ϕ metab = 0.5) and a scaling factor of f sc = 1.8 is used [1]. The metabolic component fully parallels that used for insulin secretion (n ins2,gluc = 2.5, C Hf,ins2,gluc = 7 mM; see eq. A6 later). With this selection, oxygen consumption increases about 70% when going from low glucose (3 mM) to high glucose (15 mM). A step-down function, δ, is also used to account for necrosis and cut the oxygen consumption in tissues where the oxygen concentration c oxy falls below a critical value, C cr,oxy = 0.1 µM (corresponding to p cr,oxy = 0.07 mmHg). To ensure a smooth transition, COMSOL Multiphysics' smoothed Heaviside function flc1hs, which has a continuous first derivative and no overshoot [12], is used as step-down function, δ(c oxy > C cr,oxy ) = flc1hs(c oxy -1.0×10 -4 , 0.5×10 -4 ).

Glucose consumption
Similar to oxygen consumption, glucose consumption is also assumed to follow simple Michaelis-Menten kinetics (n gluc = 1) with R max,gluc = -0.028 mol⋅m -3 ⋅s -1 and C Hf,gluc = 10 μM ( While these values are draft estimates only, they have only relatively little effect on the calculated insulin secretion. Changes in glucose concentrations due to glucose consumption by islets have only minimal influence on insulin release or cell survival because oxygen diffusion limitations in tissue or in media are far more severe than for glucose [13,14]. Even if oxygen is consumed at approximately the same rate as glucose on a molar basis and has a three-to fourfold higher diffusion coefficient (see D oxy vs. D gluc values above), this is more than offset by the differences in the concentrations available under physiological conditions. The solubility of oxygen in culture media or in tissue is much lower than that of glucose making the available oxygen concentrations much more limited (e.g., around 0.05-0.2 mM vs. 3-15 mM assuming physiologically relevant conditions) [14]. Consequently, while oxygen limitations can seriously affect islet function, glucose consumption by islets has only limited effects on the glucose levels reaching the glucose-sensing β-cells.

Insulin release
A crucial part of the model is the function describing the glucose-dependence of the insulin secretion rate, R ins . Glucose (or oxygen) is not a substrate per se for insulin production; hence, there is no direct justification for the use of Michaelis-Menten-type enzyme kinetics.
Nevertheless, the corresponding generalized form of equation A2 (Hill equation) was found to fit well the experimental results, but a Hill coefficient larger than unity (n > 1) is needed because glucose-insulin response is clearly more abrupt than the rectangular hyperbola of the Michaelis-Menten equation corresponding to n = 1 [1]. Accordingly, the main function used here to describe the glucose-insulin dynamics of the second-phase response is: with n ins2,gluc = 2.5, C Hf,ins2,gluc = 7 mM, and R max,ins2 = 3.0×10 -5 mol⋅m -3 ⋅s -1 (Table 1). These values were obtained [1] by fitting the insulin release data of human islets measured by Henquin and co-workers in staircase experiments [9]. R max,ins2 as used here corresponds to a maximum (second phase) secretion rate of ~20 pg/IEQ/min for human islets [9,15,16].
The first-phase response is incorporated into the model via a component that depends on the change (time-gradient) of glucose concentration (c t = ∂c gluc /∂t). This is non-zero only when the glucose concentration is increasing, i.e., only when c t > 0. A Hill-type sigmoid response is assumed here too: with n ins1,gluc = 2, Ct Hf,ins1,gluc = 0.03 mM⋅s -1 , and R max,ins1 = 21.0×10 -5 mol⋅m -3 ⋅s -1 (Table 1) With all these, total insulin release is obtained as the sum of first-and second-phase releases and an additional modulating function to account for the limiting effect of oxygen availability, which can become important in the core of avascular islets (especially under hypoxic conditions): An abrupt Hill-type modulating function is used as ϕ i,o (c oxy ) with n ins,oxy = 3 and C Hf,ins,oxy = 3 μM (p Hf,ins,oxy = 2 mmHg) so that insulin secretion starts becoming limited for local oxygen concentrations that are below ~6 μM (corresponding to a partial pressure of p O 2 ≈ 4 mmHg).
Finally, correct time-scaling of insulin-release could be achieved only with introduction of an extra compartment [1]; without this, insulin responses decreased too quickly compared to experimental observations (~1 min vs. ~5-10 min). Hence, insulin is assumed to be first secreted in a 'local' compartment ( Figure 1) in response to the current local glucose concentration (R ins , eq. A9) and then released from there following a first order kinetics, dc insL /dt = R ins -k insL (c insLc ins ). 'Local' insulin is modeled as an additional concentration with the regular convection model (eq. A1), but having a very low diffusivity (D insL,t = 1.0×10 -16 m 2 ⋅s -1 ). The original model, which was calibrated for human islets, used a corresponding rate constant of k insL = 0.003 s −1 , corresponding to a half-life t 1/2 of approximately 4 min [1]. Here, to fit the data obtained with murine islets, this was increased to 0.006 s −1 to have a slightly faster release; this was the only parameter modified compared to the original model.

Fluid dynamics
To incorporate media flow, these convection and diffusion models are coupled to a fluid dynamics model. The incompressible Navier-Stokes model for Newtonian flow (constant viscosity) is used for fluid dynamics to calculate the velocity field u that results from convection [2,3]: Here, ρ denotes density Computations were done with the Pardiso direct solver as linear system solver with an imposed maximum step of 0.5 s, which was needed to not miss changes in the incoming glucose concentrations that could be otherwise overstepped by the solver. Calculations are done with two spherical islets of 100 and 150 µm diameter placed in a 2D cross-section of a cylindrical tube with fluid flowing from left to right ( Figure 3 illustrates the geometry used). Islets were considered homogeneous inside; individual cells (e.g., αor β-cells) were not considered separately. Islet sizes were selected based on our analysis of the size distribution of isolated islets [17]. This confirmed that the expected value of islet diameter is 95 µm, whereas the expected value of islet volume is 1.2×10 6 µm 3 , corresponding to the volume of an islet with d = 133 µm [17]. As initial conditions, atmospheric oxygen (0.200 mol⋅m −3 ), low glucose (3 mM), and zero insulin concentrations were assumed inside the perifusion chamber. Stepwise increments in the incoming glucose concentration were implemented using the smoothed Heaviside step function at predefined time points t i , c gluc = c low + Σc step,i ·flc1hs(t -t i , τ). For FEM, COMSOL's predefined 'Extra fine' mesh size was used. In the convection and diffusion models, the following boundary conditions were used: insulation/symmetry, n⋅(-D∇c+cu) = 0, for walls, continuity for islets. For the outflow, convective flux was used for insulin, glucose, and oxygen, n⋅(-D∇c) = 0. For the inflow, inward flux was used for all components with zero for insulin (N 0 = 0), c gluc ·v in for glucose, and c oxy,in ·v in for oxygen. In the incompressible Navier-Stokes model, no slip (u = 0) was used along all surfaces corresponding to liquid-solid interfaces. For the outlet, pressure, no viscous stress with p 0 = 0 was imposed.