Pharmacokinetic Model of the Transport of Fast-Acting Insulin From the Subcutaneous and Intradermal Spaces to Blood

Abstract

Pharmacokinetic (PK) models describing the transport of insulin from the injection site to blood assist clinical decision making and are part of in silico platforms for developing and testing of insulin delivery strategies for treatment of patients with diabetes. The ability of these models to accurately describe all facets of the in vivo insulin transport is therefore critical for their application. Here, we propose a new model of fast-acting insulin analogs transport from the subcutaneous and intradermal spaces to blood that can accommodate clinically observed biphasic appearance and delayed clearance of injected insulin, 2 phenomena that are not captured by existing PK models. To develop the model we compare 9 insulin transport PK models which describe hypothetical insulin delivery pathways potentially capable of approximating biphasic appearance of exogenous insulin. The models are tested with respect to their ability to describe clinical data from 10 healthy volunteers which received 1 subcutaneous and 2 intradermal insulin injections on 3 different occasions. The optimal model, selected based on information and posterior identifiability criteria, assumes that insulin is delivered at the administrative site and is then transported to the bloodstream via 2 independent routes (1) diffusion-like process to the blood and (2) combination of diffusion-like processes followed by an additional compartment before entering the blood. This optimal model accounts for biphasic appearance and delayed clearance of exogenous insulin. It agrees better with the clinical data as compared to commonly used models and is expected to improve the in silico development and testing of insulin treatment strategies, including artificial pancreas systems.

Keywords

biphasic kinetics delayed clearance fast-acting insulin pharmacokinetic modeling

Type 1 diabetes mellitus (T1DM) is a disease of autoimmune destruction of the beta cells in the pancreas, which leads to a progressive and, ultimately, to an almost complete lack of endogenous insulin. In patients with T1DM, exogenous insulin is required to maintain blood glucose homeostasis. Accurate delivery of insulin is critical since overdosing may lead to severe and potentially fatal hypoglycemia,¹ while under-dosing could result in hyperglycemia which can be harmful both acutely and chronically.² However, accurate dosing depends on a reliable prediction of insulin appearance in, and clearance from, the blood. Such prediction is needed for the design of accurate insulin delivery profiles and is critical for the development of in silico simulation platforms for development and testing of insulin infusion algorithms.

Several routes for insulin delivery have been used clinically, including the intravenous (IV), subcutaneous (SC), intradermal (ID), interperitoneal, and, more recently, the intrapulmonary and intranasal routes.^3,4 The kinetic modeling of SC delivery has been the focus of extensive research.^5
-9 A 3-compartment model was first proposed for multiple SC insulin injections,⁹ which was later adopted for intersubject evaluation¹⁰ (see Nucci and Cobelli⁸ for a comprehensive review of the existing continuous and discrete kinetic models of SC delivery). In Wilinska et al⁵ the authors selected an optimal model over 11 candidates, in which the slower transport accounts for transition and equilibration processes between the different insulin aggregation states (monomer, dimer, hexamer).^6,7,11 Some of these kinetic models are critical components in comprehensive metabolism modeling to aid the development of insulin delivery algorithms.^12
-14 In such applications, the model parameters represent individual characteristics and, as discussed in Heinemann,¹⁵ may have significant interpatient variability due to differences in injection sites,^16,17 lipodystrophy variations,¹⁸ and/or demographic diversities.¹⁹ Accurate representation of the parameter variation has allowed for the creation of an in silico population for preclinical trials, which has been accepted by the Food and Drug Administration (FDA) as a substitute of preclinical trials.^20,21

To our knowledge, the kinetics of the ID delivery route has never been described mathematically and the goal of this study is to develop a unifying model that can accurately depict the specifics of both the ID and SC insulin delivery routes. This effort is also justified by human trials outcomes in which the ID route showed faster absorption and action of exogenous insulin as compared to the SC route.^22
-26 In addition, even though existing SC models approximate fairly well most of the in vivo insulin kinetics, certain phenomena that seem to be present in the clinical data are not represented such as the delayed spiking,^22,27 the biphasic appearance and delayed clearance.^5,28

Methods

Experimental Protocol and Clinical Data

The data were collected during a previous clinical safety and feasibility trial.²² Ten healthy male volunteers (age 31.0 ± 6.3 years, BMI 23.7 ± 1.7 kg/m², weight 75.2 ± 5.3 kg) were given a 10-unit injection of insulin analog lispro (Humalog^®, Eli Lilly, Indianapolis, IN) under euglycemic glucose clamp conditions and after an overnight fast. Each volunteer received 1 SC and 2 ID treatments on 3 separate test days with a washout period of 3-20 days. The data were collected via a very specific radioimmunoassay to the injected insulin analog with a 2% relative error of the measurement.

There are 30 data sets in total (10 SC and 20 ID profiles) and each data set consists of 17 data points. As shown in Figure 1, the insulin appearance and clearance in plasma seem to be more rapid in the ID route than in the SC route for the first 100 minutes. The 2 delivery routes also appear to have similar long-tail clearance rate. Another interesting phenomenon observed in the data is the possible biphasic appearance (2 peaks) in the bloodstream of exogenous insulin (Figure 2). The standard model may not capture accurately the observed variability in the data.

Figure 1.

Plasma insulin concentrations for SC, ID 1, and ID 2 routes on 10 subjects (mean ± SD).

Figure 2.

Biphasic phenomenon in subject 4 (SC route). Crossing symbols are the experimental data, and the red curve is the fit using a standard model as a reference.¹³

Volume of Distribution

The kinetics of insulin in plasma is determined by a single compartment model²⁹ with the assumption that all injected insulin $(\sum u (t))$ appears in the circulation and is removed at a rate of $k_{c} \min^{- 1}$ from time $t_{1}$ to $t_{2}$ , yields the equation

\int_{t_{1}}^{t_{2}} I^{'} (t) d t = - k_{c} \int_{t_{1}}^{t_{2}} I (t) d t + \frac{\sum u (t)}{V_{d}} .

Given enough time, the left-hand-side of Eq. 1 will be equal to 0. So:

V_{d} = \frac{\sum u (t)}{k_{c} \int_{t_{1}}^{t_{2}} I (t) d t} .

Model Representations

We propose for comparison 9 pharmacokinetic (PK) models of insulin appearance summarized in Table 1 ordered by level of complexity.

Table 1.

Models Representations.

#	Model structure schematics	i) Differential equations or transfer function representation(s)	Number of parameters
1		$\begin{array}{l} {\dot{x}}_{1} (t) = - (k_{d} + k_{a 1}) x_{1} (t) + u (t) \\ {\dot{x}}_{2} (t) = k_{d} x_{1} (t) - k_{a 2} x_{2} (t) \\ {\dot{x}}_{3} (t) = k_{a 1} x_{1} (t) + k_{a 2} x_{2} (t) - k_{c} x_{3} (t) \end{array}$	4
2		$\begin{array}{l} {\dot{x}}_{1} (t) = - (k_{d 1} + k_{a 1}) x_{1} (t) + u (t) \\ {\dot{x}}_{2} (t) = k_{d 1} x_{1} (t) - k_{d 2} x_{2} (t) \\ {\dot{x}}_{3} (t) = k_{d 2} x_{2} (t) - k_{a 2} x_{3} (t) \\ {\dot{x}}_{4} (t) = k_{a 1} x_{1} (t) + k_{a 2} x_{3} (t) - k_{c} x_{4} (t) \end{array}$	5
3		$\begin{array}{l} {\dot{x}}_{1} (t) = - (k_{d 1} + k_{a 1}) x_{1} (t) + u (t) \\ {\dot{x}}_{2} (t) = k_{d 1} x_{1} (t) - (k_{d 2} + k_{d 3} \frac{{(x_{2} / k_{d 4})}^{k_{d 5}}}{1 + {(x_{2} / k_{d 4})}^{k_{d 5}}}) x_{2} (t) \\ {\dot{x}}_{3} (t) = (k_{d 2} + k_{d 3} \frac{{(x_{2} / k_{d 4})}^{k_{d 5}}}{1 + {(x_{2} / k_{d 4})}^{k_{d 5}}}) x_{2} (t) - k_{a 2} x_{3} (t) \\ {\dot{x}}_{4} (t) = k_{a 1} x_{2} (t) + k_{a 2} x_{3} (t) - k_{c} x_{4} (t) \end{array}$	7 $k_{d 5} = 3.5$
3a	Same as model 3, $k_{d 5}$ is free	Same as model 3	8
4		$\begin{array}{l} {\dot{x}}_{1} (t) = - (k_{d} + k_{a 1}) x_{1} (t) + A u (t) + (1 - A) u (t - Δ) \\ {\dot{x}}_{2} (t) = k_{d} x_{1} (t) - k_{a 2} x_{2} (t) \\ {\dot{x}}_{3} (t) = k_{a 1} x_{1} (t) + k_{a 2} x_{2} (t) - k_{c} x_{3} (t) \end{array}$	6
5		Transfer function: $TF = {(\frac{b}{s + b})}^{a}$	4
6		Transfer functions: ${TF}_{1} = {(\frac{b_{1}}{s + b_{1}})}^{2} {,TF}_{2} = {(\frac{b_{2}}{s + b_{2}})}^{a_{2}}$	5
6a		Same as model 6	6
7		Transfer functions: ${TF}_{1} = {(\frac{b_{1}}{s + b_{1}})}^{2} {,TF}_{2} = {(\frac{b_{2}}{s + b_{2}})}^{a_{2}}$	6
8		Transfer functions: ${TF}_{1} = {(\frac{b_{1}}{s + b_{1}})}^{2} {,TF}_{2} = {(\frac{b_{2}}{s + b_{2}})}^{a_{2}} {,TF}_{3} = {(\frac{b_{3}}{s + b_{3}})}^{a_{3}}$	8
9		Transfer functions: ${TF}_{1} = {(\frac{b_{1}}{s + b_{1}})}^{2} {,TF}_{2} = {(\frac{b_{2}}{s + b_{2}})}^{a_{2}}$	7
10		$\begin{array}{l} {\dot{Q}}_{1 a} (t) = k \cdot u - k_{a 1} \cdot Q_{1 a} - {LD}_{a} \\ {\dot{Q}}_{1 b} (t) = (1 - k) \cdot u - k_{a 2} \cdot Q_{1 b} - {LD}_{b} \\ {\dot{Q}}_{2} (t) = k_{a 1} \cdot Q_{1 a} - k_{a 1} \cdot Q_{2} \\ {\dot{Q}}_{3} (t) = k_{a 1} \cdot Q_{2} + k_{a 2} \cdot Q_{1 b} - k_{e} \cdot Q_{3} \\ {LD}_{a} = V_{\max} \cdot Q_{1 a} / (K_{M} + Q_{1 a}) \\ {LD}_{b} = V_{\max} \cdot Q_{1 b} / (K_{M} + Q_{1 b}) \end{array}$	7

Model 1 is a standard (reference) model, currently adopted in the preclinical simulation platform approved by FDA^20,21 and for insulin treatment strategies design.^30,31 It has 3 compartments and the input is delivered at compt. 1, then it follows either (1) direct diffusion to the plasma (compt. 3) or (2) traveling through an intermediate compartment (compt. 2) before arriving at compt. 3.

Models 2, 3, and 3a are derived from model 1 with an auxiliary intermediate compartment to create a second delayed insulin spike in plasma. In model 2, all diffusions between the compartments are linear; while in model 3, a sigmoid function is used to describe the diffusion between compartments 2 and 3. In model 3, the sigmoid index is fixed $(k_{d 5} = 3.5)$ , while in model 3a this index is estimated with the other parameters.

Model 4 (delayed appearance model) imposes a direct time delay to model 1. It assumes that a fraction (A) of the input appears immediately at the administrative site (compt. 1) while the rest (1-A) arrives at comp 1 with a delay of Δ minutes.

Model 5 is described as a transfer function of the Laplace-domain variable approximating a series of diffusions from the administrative site (comp 1) to the plasma (comp 2). The parameters of the transfer function are: a (count parameter, a≥1) representing the number of diffusion compartments and b (rate parameter, min^-1) representing the diffusion rate between the compartments.

Models 6 through 9 use 2 or 3 (instead of 1, as in model 5) serial-compartmental diffusion blocks to depict various transport kinetics. In these models, insulin is delivered at the administrative site (comp 1) and is then transported through the diffusion blocks. To avoid the nonidentifiability during optimization, 1 of the diffusion rates is fixed at 0.04 min^-1, as it was shown to be nonidentifiable in preliminary analysis (model 6a, not included in analysis). The other diffusion parameters account for the distribution difference among the diffusion blocks. The count parameter $a_{1}$ of the first transfer function block is fixed at 2 to represent a fast delivery with fewer compartments. Additional intermediate compartments are incorporated in models 7 and 9, with diffusion constants k_a, k_a1, and k_a2.

Model 10 is the optimal model in Wilinska et al.⁵ The injected insulin is diverted into the slow and fast channels respectively, which have local degradations (following Michaelis-Menten relation).

Model Fitting

All 9 candidate models were fit to the 30 available data sets for each 1 of the 3 delivery profiles. As the decay of insulin during the “tail” represents a relatively slow removal process, a priori distribution is imposed on the plasma clearance parameter. Given that the insulin clearance rate (k_c) follows a log-normal distribution^32,33 and the half-life of the second phase of exogenous insulin in plasma of normal volunteers is 8.2 ± 2.2 minutes,³⁴ the normal mean (μ) and variance (σ) of the half-life are derived from the log-normal mean (m) and variance (υ) as follows:

\begin{array}{l} m = e^{μ + σ^{2} / 2} = 8.2, v = e^{2 μ + σ^{2}} (e^{σ^{2}} - 1) = {2.2}^{2} \\ μ = \log (\frac{m^{2}}{\sqrt{v + m^{2}}}) \approx 2.0694, σ = \sqrt{\log (\frac{v}{m^{2}} + 1)} \approx 0.2636 . \end{array}

Therefore, the cost function is defined as:

\begin{array}{l} f = {\sum_{i} (x_{i} - d_{i})}^{2} + {(\frac{{\hat{k}}_{c} - μ}{σ})}^{2} \\ {\hat{k}}_{c} = \ln (\frac{\ln 2}{k_{c}}) \end{array}

where x_i and d_i represent the model fit and the actual clinical data, respectively and ${\hat{k}}_{c} (μ, σ^{2}) .$ is the transformation of insulin clearance rate in plasma (k_c).

For models 5 through 9, the count parameters ( $a a$ ) are defined in stepwise loops to find the other optimal parameters. The software we are using is Matlab 2012a.

Model Selection Criteria

The Bayesian information criterion (BIC) is adopted for model selection given as

B I C = N \cdot \ln (\frac{R S S}{N}) + k \cdot \ln (N),

where N is the number of data points (excluding the initial datum at time 0), k is the degree of the optimization space, and RSS is the value of the cost function. A smaller BIC value implies a better model representation.

A 2-step procedure was implemented for model selection/rejection:

A model is rejected if its BIC values in at least 1 route are significantly larger than the corresponding BIC values of another model from the same group and the remaining BIC values are not significantly different. Priority is given to rejections based on the highest number of routes in which BIC is significantly different. A 2-tailed paired t test is used to compare the BIC values; P < .05 is considered statistically significant.

A model is rejected if 1 or more of its parameters are not identifiable. The identifiability of a parameter is determined by the Jacobian derivation method for each subject in each delivery route. A parameter will be identifiable if the value of CV is less than 100% indicating that parameter estimate accuracy is acceptable considering a posteriori identifiability.³⁵ For models 5 through 9, after finding the optimal count parameters ( $a$ ), we search the parameters space again to calculate CV values for the other parameters.

Results

After all data sets were fitted with the 9 candidate models, we found that the standard model and its derivatives (models 1 through 3) cannot reproduce biphasic phenomena due to limitations of their structure. The delayed appearance model (model 4) and the serial-compartmental diffusion models (models 5 through 9) on the other hand can approximate a biphasic kinetics. This notion is exemplified in Figures 3 and 4.

Figure 3.

Subject 4, subcutaneous route. Measured data are represented by the black crosses; all model fits are represented as continuous curves.

Figure 4.

Subject 2, intradermal route 2. Measured data are represented by the black crosses; all model fits are represented as continuous curves.

To visualize the performance of all tested models, the average residual errors of the 10 volunteers were compared in Figure 5 in 3 time intervals 0-60 minutes (immediate appearance), 60-120 minutes (late appearance), after 120 minutes (tail clearance) to test the ability to accurately capture different phases of the appearance of exogenous insulin in the bloodstream. In particular, the figure shows that models 1 through 5 fail to capture well the first 2 phases where the biphasic phenomenon occurs. The BIC values of all 9 model candidates are given in Table 2.

Figure 5.

Average absolute residual errors of the 10 subjects (3 routes; 30 data sets total) across all the model candidates for insulin appearance onset, peak and decay.

Table 2.

Average BIC Values in Model Candidates (Mean ± SD).

#	SC	ID 1	ID 2
1	72.63 ± 14.15	62.61 ± 15.61	66.86 ± 17.11
2	62.06 ± 11.39	62.08 ± 14.18	62.48 ± 15.66
3	70.60 ± 11.86	70.39 ± 15.90	68.80 ± 15.75
3a	72.47 ± 12.85	70.58 ± 13.96	71.31 ± 15.88
4	72.51 ± 14.73	70.33 ± 15.69	75.12 ± 17.35
5	58.54 ± 9.63	63.20 ± 13.82	60.71 ± 13.60
6	52.19 ± 10.42	60.84 ± 12.40	56.54 ± 16.05
6a	52.33 ± 10.75	60.25 ± 12.59	55.97 ± 15.66
7	50.54 ± 9.26	57.22 ± 12.40	55.21 ± 14.39
8	49.81 ± 9.87	54.15 ± 15.84	54.15 ± 17.42
9	51.82 ± 7.62	60.01 ± 12.22	59.15 ± 14.35
10	68.09 ± 12.48	66.67 ± 16.02	65.84 ± 16.80

Model Selection

Models 3, 3a, 4, and 10 are rejected first since their BIC values were significantly larger in all routes than the BIC values of model 7. Moreover, model 3a (with the free index) has significant greater BIC values than those of model 3 (with the fixed index) in Route SC and ID2. Models 1, 2, 5, and 9 are rejected since 2 of their BIC values are significantly larger than the corresponding BIC values of model 7 and the remaining BIC values are not significantly different. Analogously, we reject models 6 and 6a which are outperformed by model 7 in the ID 1 route respectively. Among the remaining 2 model candidates, models 7 and 8, model 8 are rejected since some of their parameters are nonidentifiable (CV > 200%). The model selection/rejection process is recapitulated in Table 3. The parameters of model 7 are identifiable (CV < 100%) shown in Table 4 (mean ± SD).

Table 3.

Model Selection Process.

Step	Comparison	Primary criterion (BIC)			Secondary criterion	Rejected model
Step	Comparison	SC	ID1	ID2	(Identifiability)	Rejected model
1	Model 7 vs model 3	P < .0001	P = .0025	P = .0005	N/A	Model 3
2	Model 7 vs model 3a	P < .0001	P = .0004	P = .0001	N/A	Model 3a
3	Model 7 vs model 4	P = .0003	P = .0028	P = .0001	N/A	Model 4
4	Model 7 vs model 10	P = .0010	P = .0024	P = .0001	N/A	Model 10
5	Model 7 vs model 1	P = .0003	NS	P = .0029	N/A	Model 1
6	Model 7 vs model 2	P = .0047	NS	P = .022	N/A	Model 2
7	Model 7 vs model 5	P = .01	NS	P = .025	N/A	Model 5
8	Model 7 vs model 9	NS	P = .0003	P < .0001	N/A	Model 9
9	Model 7 vs model 6	NS	P = .033	NS	N/A	Model 6
10	Model 7 vs model 6a	NS	P = .011	NS	N/A	Model 6a
11	Model 7 vs model 8	NS	NS	NS	Nonidentifiable	Model 8
12	Model 7	Not applicable			Identifiable	—

Table 4.

Average Optimal Parameter Values of Model 7 (Mean ± SD) and CV Values.

Route	a₂	k_d	b₁	b₂	k_a	$\hat{k}$ _c
SC	15 ± 14 (N/A)	0.146 ± 0.130 (9.77%)	0.088 ± 0.068 (10.43%)	0.227 ± 0.162 (6.64%)	0.015 ± 0.012 (10.57%)	2.095 ± 0.246 (10.29%)
ID1	17 ± 15 (N/A)	0.160 ± 0.094 (15.52%)	0.217 ± 0.119 (12.41%)	0.574 ± 0.514 (10.38%)	0.019 ± 0.028 (6.01%)	2.183 ± 0.461 (8.82%)
ID2	9 ± 6 (N/A)	0.138 ± 0.098 (10.39%)	0.307 ± 0.372 (12.09%)	0.327 ± 0.238 (6.75%)	0.015 ± 0.010 (4.78%)	2.500 ± 0.504 (5.92%)

In summary, model 7 is selected as the optimal choice among all candidates. Its average parameters are given in Table 4. Model 7 recreated similar insulin appearance shown in Figure 1, which demonstrates the excellent capability of the optimal model to track the PK specifics of the experimentally observed exogenous insulin appearance and clearance. The residuals are demonstrated in Figure 6 as well.

Figure 6.

Average residual errors across all the model candidates for each insulin concentration measurement across all routes and subjects.

Discussion

This work uses data from a previous clinical safety and feasibility trial²² to develop and verify a new unified PK model of a fast-acting insulin analog transport from the subcutaneous and intradermal spaces to the circulation. This model addresses a need for an accurate prediction of insulin transport that accommodates a biphasic appearance and delayed clearance of exogenous insulin required for the development of in silico simulation platforms and insulin delivery algorithms.

To develop and validate the new model we have adopted the following approach. Nine models with different complexities are proposed for assessment including 3 standard models (models 1, 2, 3, and 3a), a delayed appearance model (model 4), and 5 serial-compartmental diffusion models (models 5 through 9). The performance of all models is estimated based on their ability to describe the available clinical data. The individual fits optimize a cost function modified to penalize deviations of the elimination rate k_c from a value corresponding to a literature derived half-life of circulating insulin. This was done to guarantee that the models would capture accurately the final phase of insulin appearance in the bloodstream (the “tail”: 120-360 min) which is needed for safe model application (below).

In terms of model performance, models 1, 2, 3, and 3a had larger errors in the first 2 hours after the infusion (Figure 3) corresponding failure to capture correctly the tendency in the data of an insulin spike to rapidly appear within an hour after the injection followed by a second delayed spike during the second hour after the injection. These large errors could indicate that the structure of the standard model and its derivatives is not “rich” enough to reproduce the biphasic kinetics of insulin which appears present in some of the data sets. The underlying cause for the differences in the PK parameters between the SC and ID delivery routes has not been completely established, but may be due to differences in the local fat tissue structure¹⁹ or to differences in the vasculature of the infusion site with a tendency for the venous capillaries to mediate the insulin transport from the SC tissue and for the lymphatic vessels to facilitate the appearance after an ID injection.²²

Model 4 adopts a delayed release and imposes immediate distribution of injected insulin into 2 separate channels, but still may not represent the first and the second peaks (Figures 3). Similarly, model 5 (single serial-compartmental diffusions) is not able to depict well the biphasic phenomenon either.

Models 6, 6a, 7, 8, and 9 adopt different strategies to compensate for the limitation of the single serial-compartmental diffusions. Models 6 and 6a have dual serial-compartmental diffusions. Model 7 adds an intermediate diffusion compartment, model 8 adds a third diffusion compartment, and model 9 adds 2 intermediate diffusion compartments before plasma. These 4 models improved the fitting showing that their structures provides a better approximation of the physiology of in vivo exogenous insulin delivery. Among these models, model 7 is chosen as the best representation for an insulin PK transport model.

To summarize, models 1, 2, and 3 cannot reproduce a dual phase insulin appearance. Models 4 and 5 have larger errors over the population as compared to models 6, 6a, 7, 8, and 9 which can reproduce a biphasic appearance and a slower final clearance phase. The values of BIC will confirm that increasing the model complexity and the degree of optimization space does not interfere with the ability of the models to describe the data. to avoid the overfitting of the limited points in the current data set, we did not investigate models more complex than models 7 to 9. The physiological relevance of model 7 is also supported by the identifiability of its parameters.

In 1 application of the current modeling efforts, we tested whether there exists a relationship between a prolonged clearance and the delayed appearance of the insulin spike(s). In particular, model 7 was used to analyze all 10 volunteers in the 3 delivery routes separately. The last insulin peak time were recorded and compared to the corresponding fitted half-lives in plasma. We found that there was no significant association (r < .5, P > .5) between the onset time of the last insulin peak and the half-life in plasma. However, this may not be sufficient to suggest the independence between the insulin transport in tissue and its removal from plasma since a weighting of half-life is defined in the cost function, which may need further study.

Our results demonstrate that models capable of recreating a biphasic appearance of exogenous insulin agree better with the clinical data. In particular, the selected optimal unified (SC & ID) PK model significantly outperforms existing standard models, which cannot recreate biphasic appearance. This outcome supports the concept that the biphasic phenomenon exists and suggests a putative system-level physiological mechanism that could mediate it. This mechanism is embodied in the optimal chosen model (model 7)—a serial-compartmental diffusion model that assumes that insulin is delivered at the administrative site and is further transported to the circulation through 2 separate diffusion routes.

Although biphasic appearance seems present in the individual data from many volunteers in all administration routes, the overall population averages shown in Figure 1 does not have this property due to significant intersubject variability, such as widely distributed time and/or magnitude of the onset of the second phase appearance. The data in Figures 2A and 2B show different insulin profiles (2 distinct, widely separated phases and comparable dual-peak levels vs 2 peaks close to each other and at different levels). This large intersubject variability favors future efforts to concentrate on creating individual distributions for use in in silico simulation platforms rather than using averaged population parameters.

A limitation of the current study is that 10 subjects is not enough to consider the model valid across a broader population. Moreover, the data were obtained from normal volunteers. Some physiological characteristics approximated by the model may be different in patients with T1DM and it may be necessary to create parameter distributions specific for T1DM patients before the new model can be fully incorporated into in silico pretrial simulation platforms (such as the UVa/Padova T1DM simulator).

Conclusion

Existing standard compartmental linear models of insulin transport form the injection site to the circulation can reproduce only a single phase mode of insulin appearance. Clinical data in healthy volunteers however suggest a mode of biphasic appearance after SC or ID delivery. This work proposes an insulin transport model based on the comparison of 9 model candidates with increasing complexities deliberately introduced to depict the biphasic phenomenon. An information criteria-based analysis shows that a model with 2 serial-compartmental diffusions and 1 intermediate compartment provides the best description of the clinical data in both the SC and ID delivery routes. The selected optimal model can be further incorporated into control algorithms and simulating platforms, or can be used to assist the development of more accurate profiles of exogenous insulin appearance with an immediate clinical impact.

Footnotes

Acknowledgements

The authors would like to thank BD for providing the insulin data, Stephen Patek and Edward Ortiz for physiology discussions, and Stacey Anderson for discussion on clamp procedures.

Abbreviations

BIC, Bayesian information criterion; FDA, Food and Drug Administration; ID, intradermal; IV, intravenous; PK, pharmacokinetics; SC, subcutaneous; T1DM, type 1 diabetes mellitus.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported in part by the National Institute of Health under grants R01 DK51562, R01 DK082805, and JDRF 17-2011-273

References

Cryer

. Hypoglycemia in diabetes: pathophysiological mechanisms and diurnal variation. Progress Brain Res. 2006;153:361-365.

Klein

. Hyperglycemia and microvascular and macrovascular disease in diabetes. Diabetes Care. 1995;18(2):258-268.

Owens

Zinman

Bolli

. Alternative routes of insulin delivery. Diabetic Med. 2003;20(11):886-898.

Smith

Sargent

Wilson

Price

. Subcutaneous or intramuscular insulin injections. Arch Dis Childhood. 1991;66(7):879-882.

Wilinska

Chassin

Schaller

Schaupp

Pieber

Hovorka

. Insulin kinetics in type-1 diabetes: continuous and bolus delivery of rapid acting insulin. IEEE Trans Biomed Eng. 2005;52(1):3-12.

Tarín

Teufel

Picó

Bondia

Pfleiderer

. Comprehensive pharmacokinetic model of insulin Glargine and other insulin formulations. IEEE Trans Biomed Eng. 2005;52(12):1994-2005.

Kuang

. Systemically modeling the dynamics of plasma insulin in subcutaneous injection of insulin analogues for type 1 diabetes. Math Biosci Eng. 2009;6(1):41-58.

Nucci

Cobelli

. Models of subcutaneous insulin kinetics. A critical review. Computer Methods Programs Biomed. 2000;62(3):249-257.

Puckett

Lightfoot

. A model for multiple subcutaneous insulin injections developed from individual diabetic patient data. Am J Physiol Endocrinol Metab. 1995;(269):E1115-E1124.

10.

Clausen

De Gaetano

Vølund

. Within-patient variation of the pharmacokinetics of subcutaneously injected biphasic insulin aspart as assessed by compartmental modeling. Diabetologia. 2006;49(9):2030-2038.

11.

Kang

Brange

Burch

Vølund

Owens

. Absorption kinetics and action profiles of subcutaneously administered insulin analogues (AspB9GluB27, AspB10, AspB28) in healthy subjects. Diabetes Care. 1991;14(11):1057-1065.

12.

Cobelli

Dalla Man

Sparacino

. Diabetes: models, signals, and control. IEEE Rev Biomed Eng. 2009;2:54-96.

13.

Dalla Man

Raimondo

Rizza

Cobelli

. GIM, simulation software of meal glucose-insulin model. J Diabetes Sci Technol. 2007;1(3):323-330.

14.

El-Khatib

Russell

Nathan

Sutherlin

Damiano

. A bihormonal closed-loop artificial pancreas for type 1 diabetes. Sci Translational Med. 2010;2(27):27ra27.

15.

Heinemann

. Variability of insulin absorption and insulin action. Diabetes Technol Ther. 2002;4(5):673-682.

16.

Ter Braak

Woodworth

Bianchi

. Injection site effects on the pharmacokinetics and glucodynamics of insulin lispro and regular insulin. Diabetes Care. 1996;19(12):1437-1440.

17.

Berger

Cüppers

Hegner

Jörgens

Berchtold

. Absorption kinetics and biologic effects of subcutaneously injected insulin preparations. Diabetes Care. 1982;5(2):77-91.

18.

Kølendorf

Bojsen

Deckert

. Clinical factors influencing the absorption of 125I-NPH insulin in diabetic patients. Hormone Metab Res. 1983;15(6):274-278.

19.

Sindelka

Heinemann

Berger

Frenck

Chantelau

. Effect of insulin concentration, subcutaneous fat thickness and skin temperature on subcutaneous insulin absorption in healthy subjects. Diabetologia. 1994;37(4):377-380.

20.

Kovatchev

Breton

Dalla Man

Cobelli

. In silico preclinical trials: a proof of concept in closed-loop control of type 1 diabetes. J Diabetes Sci Technol. 2009;3(1):44-55.

21.

Dalla Man

Rizza

Cobelli

. Meal simulation model of the glucose-insulin system. IEEE Trans Biomed Eng. 2007;54(10):1740-1749.

22.

Pettis

Ginsberg

Hirsch

. Intradermal microneedle delivery of insulin lispro achieves faster absorption and insulin action than subcutaneous injection. Diabetes Technol Ther. 2011;13(4):435-442.

23.

Pettis

Hirsch

Kapitza

. Microneedle-based intradermal versus subcutaneous administration of regular human insulin or insulin lispro: pharmacokinetics and postprandial glycemic excursions in patients with type 1 diabetes. Diabetes Technol Ther. 2011;13(4):443-450.

24.

Norman

Brown

Raviele

Prausnitz

Felner

. Faster pharmacokinetics and increased patient acceptance of intradermal insulin delivery using a single hollow microneedle in children and adolescents with type 1 diabetes. Pediatr Diabetes. 2013;14(6):459-465.

25.

Gupta

Felner

Prausnitz

. Rapid pharmacokinetics of intradermal insulin administered using microneedles in type 1 diabetes subjects. Diabetes Technol Ther. 2011;13(4):451-456.

26.

McVey

Hirsch

Sutter

. Pharmacokinetics and postprandial glycemic excursions following insulin lispro delivered by intradermal microneedle or subcutaneous infusion. J Diabetes Sci Technol. 2012;6(4):743-754.

27.

Steiner

Hompesch

Pohl

. A novel insulin formulation with a more rapid onset of action. Diabetologia. 2008;51(9):1602-1606.

28.

Breton

Chan

Keith

. Pharmacodynamic and pharmacokinetic analysis of micro-needle based intradermal insulin delivery using a biphasic model. Paper presented at: International Conference on Advanced Technologies & Therapies for Diabetes; 2010; Basel, Switzerland.

29.

Dhillon

Kostrzewski

. Clinical Pharmacokinetics. London, UK: Pharmaceutical Press, 2006.

30.

Abu-Rmileh

Garcia-Gabin

. A gain-scheduling model predictive controller for blood glucose control in type 1 diabetes. IEEE Trans Biomed Eng. 2010;57(10):2478-2484.

31.

Wang

Dassau

Doyle

. Closed-loop control of artificial pancreatic β-cell in type 1 diabetes mellitus using model predictive iterative learning control. IEEE Trans Biomed Eng. 2010;57(2):211-219.

32.

van Haeften

Pimenta

Mitrakou

. Disturbances in β-cell function in impaired fasting glycemia. Diabetes. 2002;51(suppl 1):S265-270.

33.

Magni

Raimondo

Bossi

. Model predictive control of type 1 diabetes: an in silico trial. J Diabetes Sci Technol. 2007;1(6):804-812.

34.

Dreyer

Matthaei

Kühnau

Rüdiger

. Prolonged plasma half-life of insulin in patients with a genetic defect of high affinity binding sites. Hormone Metab Res. 1986;18(4):247-249.

35.

Cobelli

Foster

Toffolo

. Tracer Kinetics in Biomedical Research: From Data to Model. New York, NY: Kluwer; 2000.