Modelling adrenaline secretion during counterregulatory response in Type 1 Diabetes for improved hypoglycaemia prediction

9th IFAC Symposium on Biological and Medical Systems 9th IFAC Symposium on Biological and Medicalonline Systems Available at www.sciencedirect.com Aug. 31 - Sept. 2, 2015. Berlin, Germany Aug. 31 - Sept. 2, 2015. Berlin, Germany

ScienceDirect IFAC-PapersOnLine 48-20 (2015) 213–218

Modelling Modelling adrenaline adrenaline secretion secretion during during counterregulatory counterregulatory response response in in Type Type 11 Diabetes Diabetes for for improved improved hypoglycaemia hypoglycaemia prediction prediction Vanessa Vanessa Moscardó*, Moscardó*, Paolo Paolo Rossetti Rossetti **, **, F. F. Javier Javier Ampudia-Blasco***, Ampudia-Blasco***, Jorge Jorge Bondia*† Bondia*†  

*Universitat *Universitat Politècnica Politècnica de de València, València, Valencia, Valencia, Spain Spain (( e-mail: [email protected], [email protected]). e-mail: [email protected], [email protected]). **Hospital **Hospital Sant Sant Francesc Francesc de de Borja, Borja, Gandia, Gandia, Spain, Spain, (e-mail: (e-mail: [email protected]) [email protected]) *** Hospital Clinico Universitario de Valencia,Valencia, *** Hospital Clinico Universitario de Valencia,Valencia, Spain, Spain, (e-mail: (e-mail: [email protected]) [email protected]) †Corresponding †Corresponding author author Abstract: Abstract: Counterregulatory Counterregulatory hormones hormones play play an an important important role role in in recovery recovery from from hypoglycaemia. hypoglycaemia. For For this this reason, the study of their secretion under hypoglycaemic conditions may be relevant for the prediction reason, the study of their secretion under hypoglycaemic conditions may be relevant for the prediction of of recovery recovery from from low low glucose glucose episodes. episodes. In In Type Type 11 Diabetes Diabetes Mellitus Mellitus patients, patients, adrenaline adrenaline is is the the first first line line of of counterregulatory counterregulatory hormone hormone response response to to hypoglycaemia hypoglycaemia due due to to the the early early impairment impairment of of glucagon glucagon secretion. An adrenaline secretion model to understand its secretion during hypoglycaemia is secretion. An adrenaline secretion model to understand its secretion during hypoglycaemia is proposed proposed in in this this paper, paper, with with focus focus on on the the analysis analysis of of inter-patient inter-patient variability. variability. Data Data from from hypoglycaemic hypoglycaemic clamps clamps were were used used for for the the model model development. development. Individual Individual and and population population parameter parameter values values estimation estimation was was carried carried out out by means of global optimization and Markov Chain Monte Carlo techniques, respectively. by means of global optimization and Markov Chain Monte Carlo techniques, respectively. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved.

Keywords: Keywords: Modelling, Modelling, Biomedical Biomedical systems, systems, Type Type 11 Diabetes, Diabetes, Hypoglycaemia, Hypoglycaemia, Counterregulation Counterregulation  

1. 1. INTRODUCTION INTRODUCTION The The risk risk of of hypoglycaemia hypoglycaemia is is still still aa limiting limiting factor factor of of achieving near-normoglycaemia in Type 1 Diabetes achieving near-normoglycaemia in Type 1 Diabetes Mellitus Mellitus patients. However, counterregulatory response to patients. However, counterregulatory response to hypoglycaemia is often neglected in simulation studies for hypoglycaemia is often neglected in simulation studies for control control strategies strategies design design or or predictive predictive algorithms algorithms for for hypoglycaemia prevention. This may lead to hypoglycaemia prevention. This may lead to underestimated underestimated hypoglycaemia hypoglycaemia recovery recovery and and unwanted unwanted glycaemic glycaemic rebounds. rebounds. During During hypoglycaemia, hypoglycaemia, while while the the glucose glucose concentration concentration is is decreasing, central and peripheral glucose decreasing, central and peripheral glucose sensors sensors activate activate mechanisms mechanisms to to elicit elicit the the neuroendocrine, neuroendocrine, autonomic autonomic and and behavioural response (the so-called counterregulatory behavioural response (the so-called counterregulatory response) response) (Cryer, (Cryer, 2001). 2001). In In healthy healthy humans, humans, the the initial initial response to prevent a decline in blood glucose concentration response to prevent a decline in blood glucose concentration is is aa reduction reduction in in endogenous endogenous insulin insulin secretion, secretion, which which begins begins while plasma glucose concentration is still while plasma glucose concentration is still in in the the normoglycaemic normoglycaemic range, range, at at approximately approximately 80 80 mg/dL. mg/dL. Glucagon Glucagon and and adrenaline adrenaline are are secreted secreted as as glucose glucose levels levels fall fall slightly below the normoglycaemic range, at approximately slightly below the normoglycaemic range, at approximately 68 68 mg/dL mg/dL (Schwartz (Schwartz et et al, al, 1987; 1987; Tesfaye Tesfaye and and Seaquist, Seaquist, 2010). 2010). Additionally, there is an activation of the autonomic Additionally, there is an activation of the autonomic nervous nervous system system that that increases increases adrenaline adrenaline in in the the circulation. circulation. Beyond Beyond the reduction of insulin, glucagon plays the reduction of insulin, glucagon plays the the primary primary role role in in the correction of hypoglycaemia while adrenaline the correction of hypoglycaemia while adrenaline has has aa secondary secondary role, role, as as shown shown by by experiments experiments in in which which recovery recovery from acute hypoglycaemia was studied in healthy from acute hypoglycaemia was studied in healthy subjects subjects (Cryer, (Cryer, 2002; 2002; Rizza Rizza et et al, al, 1979). 1979). Glucagon Glucagon is is secreted secreted by by

pancreatic pancreatic alpha-cells alpha-cells and and acutely acutely raises raises plasma plasma glucose glucose concentration by stimulating hepatic glucose concentration by stimulating hepatic glucose production production via via glycogenolysis glycogenolysis and and gluconeogenesis. gluconeogenesis. Adrenaline Adrenaline acts acts on on alpha alpha and and beta beta adrenergic adrenergic receptors receptors at at multiple multiple end end organs organs to to effect effect aa more more sustained sustained increase increase in in plasma plasma glucose glucose concentration: concentration: adrenaline adrenaline increases increases hepatic hepatic glycogenolysis glycogenolysis and gluconeogenesis; reduces insulin and gluconeogenesis; reduces insulin secretion secretion while while increasing increasing glucagon glucagon release release from from the the pancreatic pancreatic islets; islets; reduces reduces glucose glucose uptake uptake and and utilization utilization and and increases increases glycolysis by muscle; and, in addition, increases glycolysis by muscle; and, in addition, increases lipolysis lipolysis in in adipose adipose tissue tissue (Cryer (Cryer et et al, al, 2003). 2003). Increased Increased secretion secretion of of growth growth hormone hormone and and cortisol cortisol also also contributes contributes to to the the response response to hypoglycaemia. However, these two hormones to hypoglycaemia. However, these two hormones do do not not have have an immediate role in the recovery from hypoglycaemia an immediate role in the recovery from hypoglycaemia (Cryer (Cryer et et al, al, 1987), 1987), as as they they are are mainly mainly involved involved in in the the reaction to prolonged hypoglycaemia (Cryer and Gerich, reaction to prolonged hypoglycaemia (Cryer and Gerich, 1985; 1985; Cryer, Cryer, 2001). 2001). In In Type Type 11 Diabetes, Diabetes, counterregulatory counterregulatory response response to to hypoglycaemia is impaired. Firstly, due to exogenous hypoglycaemia is impaired. Firstly, due to exogenous insulin insulin therapy, therapy, patients patients with with Type Type 11 Diabetes Diabetes cannot cannot reduce reduce systemic insulin levels as blood glucose concentrations systemic insulin levels as blood glucose concentrations begin begin to to decline, decline, unless unless an an artificial artificial pancreas pancreas or or sensor-augmented sensor-augmented pump pump is is used used (with (with the the inherent inherent limitations limitations of of subcutaneous subcutaneous insulin delivery). Thus, subjects with Type insulin delivery). Thus, subjects with Type 11 Diabetes Diabetes lack lack the the first first line line of of defence defence against against hypoglycaemia. hypoglycaemia. Secondly, Secondly, glucagon glucagon secretion secretion in in response response to to hypoglycaemia hypoglycaemia is is lost lost soon soon after the onset of the disease. Thirdly, the response after the onset of the disease. Thirdly, the response of of adrenaline adrenaline to to aa given given level level of of hypoglycaemia hypoglycaemia is is blunted blunted and and

2405-8963 © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Copyright © 2015 213Control. Peer review underIFAC responsibility of International Federation of Automatic Copyright © 2015 IFAC 213 10.1016/j.ifacol.2015.10.141

9th IFAC BMS 214 Vanessa Moscardó et al. / IFAC-PapersOnLine 48-20 (2015) 213–218 Aug. 31 - Sept. 2, 2015. Berlin, Germany

the glycaemic threshold for its secretion is shifted to lower plasma glucose concentration (Amiel et al, 1988; Cryer et al, 2003), together with reduced autonomic symptoms. Nevertheless, due to the absence of endogenous insulin and glucagon response, adrenaline becomes the main actor of the counterregulatory response to hypoglycaemia in Type 1 Diabetes.

a biphasic nature: the first phase is a quick response against hypoglycaemia that starts when glycaemia is lower than a given threshold; the late phase is associated to the recovery presenting different dynamics. Adrenaline concentration remains constant at basal values until counterregulatory response begins.

To date, few works have partially included counterregulatory response in prediction or simulation models for a more accurate description of hypoglycaemia episodes and timecourse recovery from them. A main difficulty is the quantification of counterregulatory plasma hormonal concentrations. Kovatchev et al (1999), presented a mathematical model of insulin-glucose dynamics that included estimates for the onset and rate of counterregulatory responses. They used plasma adrenaline concentrations to validate model results, as they showed high correlation with counterregulation rates. More recently, Dalla Man et al (2014) have included in the latest version of the UVA/PADOVA Type 1 Diabetes Simulator an improved description of glucose kinetics in hypoglycaemia. They have incorporated glucagon secretion and action, as well as a paradoxical increment of glucose utilization during hypoglycaemia.

Fig. 1. Example of adrenaline response to hypoglycaemia. A parallel-input compartmental model is proposed to describe this behaviour, with equations:

Q1 (t )   k a1·Q1 (t )  1·u (t ) Q (t )  k (Q (t )  Q (t ))

In this manuscript a mathematical model of adrenaline secretion during hypoglycaemia is presented, as a first step towards modelling of the counterregulatory response in Type 1 Diabetes. The paper is organised as follows: Section 2 presents the structure of proposed mathematical model, the algorithms used to estimate parameters and the methods to evaluate the model; Section 3 shows the results of parameter estimation and goodness of fit; conclusions are drawn in Section 4.

2

a1

1

2

Q 3 (t )  ka1·Q2 (t )  ke ·Q3 (t )  ka 2 ·Q4 (t ) Q (t )   k ·Q (t )   ·u (t ) 4

a2

4

2

0 G (t )  Gth  u (t )  G (t )  Gth G (t )  Gth

(1) (2) (3) (4)

(5)

2. METHODS

A(t )  Q3 (t ) V  Abasal

(6)

2.1 Dataset

Q1 (0)  Q2 (0)  Q3 (0)  Q4 (0)  0;

Data from a eu-hypoglycaemic clamp study was used. Fourteen subjects with Type 1 Diabetes mellitus were enrolled in the study performed at the Clinic University Hospital of Valencia, Spain. Each individual participated in two eu-hypoglycaemic clamp studies with different levels of insulinemia (0.3 mU/Kg/min vs. 1 mU/Kg/min). In an initial phase glucose was normalized to 90 mg/dL by using a variable i.v. insulin infusion. Then a hypoglycaemic plateau at 50 mg/dL was induced for 45 minutes with previous and subsequent phases of euglycaemia. Total duration of the study was 8.5 hours. Plasma glucose was measured every five minutes (YSI 2300, YSI Incorporated Life Sciences, Yellow Springs, Ohio, USA). Plasma adrenaline concentration was measured every 30 minutes due to blood sampling limitations.

where A(t) is plasma adrenaline concentration (ng/L); Abasal is basal adrenaline concentration (ng/L); V is the distribution volume of adrenaline (L); G(t) is plasma glucose concentration (mg/dL); Gth is the glucose concentration value that activates the counterregulatory response of adrenaline; u(t) is the model input corresponding to glucose deviation from the activation threshold Gth..(glucose does not affect adrenaline secretion when above this threshold); Q3(t) is the measurement compartment of adrenaline mass; rest of compartments, Q1(t), Q2(t) and Q4(t), define the dynamics of each secretion phase (second and first order, respectively); β1 and β2 represent the gain of physiological response; ka1 and ka2 are transfer rate constants between compartments; and ke is adrenaline rate of disappearance.

2.2 Mathematical model

Figure 2 depicts the model structure corresponding to equations (1)-(6).

Figure 1 represents a typical response for adrenaline in our dataset after data interpolation. Adrenaline secretion presents 214

9th IFAC BMS Aug. 31 - Sept. 2, 2015. Berlin, Germany Vanessa Moscardó et al. / IFAC-PapersOnLine 48-20 (2015) 213–218

215

greater than 1 when glycaemia was lower than 70mg/dL. Gth was defined as the value of blood glucose in the time instant of change. The distribution volume of adrenaline, V, was fixed to 20L, taking as a reference studies in the scientific literature (Dejgaard et al, 1989) for an average patient. The global optimization algorithm CMA-ES (Auger and Hansen, 2005; Hansen 2005; Hansen 2006) was then used for the estimation of the parameters θ. A least squares error function was considered as the cost function to be minimized:

ˆ  argmin

Fig. 2. Compartmental model of adrenaline secretion.



i

2

i

(8)

i 1

where Aˆi ( , t ) is the predicted adrenaline value at instant t, Ai (t ) is the measured value and N is the number of data

2.2.1 Model identifiability. Consider the unknown parameter vector θ= [β1, β2, ka1, ka2, ke]T. From (1)-(6), adrenaline concentration at time ti, A(ti), is predicted by a function of the unknown parameter vector θ and the unknown input threshold Gth, i.e., A(ti)=h(ti, θ, Gth).

points. ˆ and  are the estimated parameters vector and the parameter vector respectively. A fourth-order Runge-Kutta method (ode45 in Matlab) was used for model simulation during the optimization process. By fitting individual data sets, the structure of the adrenaline secretion model was evaluated. All calculations were performed in MATLAB R2012a.

A priori (structural) identifiability deals with the uniqueness of the solution with respect to model parameters θ in the whole complex space under the ideal conditions of error-free model structure and noise-free data. Thus, to prove structural identifiability of a parameter  it is equivalent to prove if equal values of the model output A(t) imply only one (global identifiable), a finite (locally identifiable) or an infinitive (not identifiable) number of solutions in C for the corresponding parameter vector  (Carson and Cobelli, 2001).

2.3.2 Population parameter estimation Once the validity of the model structure was tested, population parameter values were obtained by means of Bayesian inference (Wilkinson, 2011). The main purpose of Bayesian inference is to determine the posterior distribution of individual and population parameters. For the provided model, this cannot be achieved by direct (analytical) computations. Instead, sampling techniques, such the Markov chain Monte Carlo (MCMC) method, have to be used. These sampling techniques provide a large sample of the posterior distributions normally running into thousands to tens of thousands of samples. The posterior distributions were summarized by the median. Computations were performed with the public domain software WinBUGS (Lunn, 2000).

There are many methods for testing models for structural identifiability. For linear systems, a common method is the analysis of the number of solutions of the equation:

y ( s, ˆ)  y ( s,  )  0

N

 ( Aˆ ( , t )  A (t ))

s, u ( s )

(7)

where y(s,.) is the system output Laplace transform (i.e. adrenaline in our case),  the (unknown) actual values of the model parameters, ˆ the fictitious model parameter values leading to the same output and u(s) the system input. The above equation was solved by comparison of transfer function coefficients leading to a set of algebraic equations (Walter and Pronzato, 1997).

Subjects were assumed to be representatives from the same homogenous population. In particular, we assumed that the individual parameters µj are drawn from a multivariate lognormal distribution:

 j  [log( 1 ), log(  2 ), log( k a1 ), log( k a 2 ), log( ke )] ~ N (  , ) 2.3 Parameter estimation (9) Parameter estimation was divided in two parts: firstly, an individual estimation and then, a population approach.

where µ is an unknown population mean vector, and Σ is an unknown covariance matrix, with all off-diagonal values equal to 0. In addition, for Bayesian analysis, specification of prior distributions of population parameters is required. It was assumed that mean of each parameter (µk) follows a normal distribution whose mean and standard deviation is based on the individual distribution of each parameter, obtained from previous results with real data; and covariance was described by:

2.3.1 Individual parameter estimation In order to reduce the number of dimensions in the global optimisation problem, prior to the identification of the model parameters vector θ, an identification of the glucose threshold (Gth) was carried out for each subject by detecting a significant change in the slope of the adrenaline temporal signal. Criterion to define a relevant change was a slope

 kk ~ Gamma(0.001,0.001) 215

k= 1,…5

(10)

9th IFAC BMS 216 Vanessa Moscardó et al. / IFAC-PapersOnLine 48-20 (2015) 213–218 Aug. 31 - Sept. 2, 2015. Berlin, Germany

Evaluating the equality (7), there is only one set of parameters θ that can lead to a given output. This means that model is structurally identifiable.

2.4 Analysis of Results The accuracy of the data fit was tested by the coefficient of determination (R2), usually interpreted as the percentage of the total variation of the dependent variable around its mean that is explained by the fitted model, and by the normalized root mean square error (NRMSE) where lower values indicate less residual variance. Therefore, the goodness of fit of the model was evaluated by the closeness of its coefficient of determination to 100% and its NRMSE to 0. The coefficient of variation of the RMSE, CV(RMSE), was also considered to observe its variability: R2  1 

VARres VARtot

1  ˆ1 ;  2  ˆ2 ; ka1  kˆa1; ka 2  kˆa 2 ; ke  kˆe 3.2 Individual parameter estimation

The glucose threshold values shown in Table 1 were obtained. In tables, numbers in brackets indicate the standard deviation (sd) of estimated values. The average value was approximately 60mg/dL (60.58±6.57 mg/dL), consistent with the findings of other studies (Schwartz et al, 1987; Tesfaye and Seaquist, 2010). Variability of this value could be caused by extra-patient variability (for example variability from noise in the data), by physiological inter-individual variability or likely by previous episodes of hypoglycaemia affecting both the degree and the threshold of the adrenergic response to hypoglycaemia (Beall and McCrimmon, 2012, de Galan et al, 2003). The identified values were consistent as demonstrated by a successful posterior identification of the model parameters θ.

(11)

N

 ( Aˆ (t )  A(t )) i

RMSE 

2

(12)

i

i 1

N RMSE NRMSE  Amax  Amin

(13)

RMSE A

(14)

CV ( RMSE ) 

(21)

Table 1. Individual estimation of Gth

where VARres and VARtot is the variance of the model residual and the observations respectively; Aˆ (ti ) is predicted value at instant ti; A(ti ) is measured value; and N is the number of predictions. 3. RESULTS Adrenaline measurements were available from 24 clamp studies. Before starting with identification, data from three clamp studies were excluded because of inaccurate adrenaline measurements (noise greater than response to hypoglycaemia).

#

Gth (mg·dL-1)

#

Gth (mg·dL-1)

#

Gth (mg·dL-1)

1

52.61 (2.80)

8

61.95(1.36)

15

58.86(1.02)

2

58.65(3.56)

9

71.27(1.02)

16

60.49(2.46)

3

52.46(3.84)

10

65.53(3.26)

17

63.74(1.89)

4

50.78(2.19)

11

63.39(2.97)

18

63.55(1.41)

5

74.53(1.55)

12

60.15(2.72)

19

49.19(3.32)

6

65.28(2.60)

13

52.87(1.40)

20

62.29(2.34)

7

67.31(1.03)

14

59.28(1.79)

21

57.95(2.20)

3.1 Identifiability Considering the state-space of the model (1)-(6):  k a1 k A   a1  0   0



C 0 0

0  k a1 k a1 0 1 V

0

0 0  ke 0

 1  0  0  0  B   2  ka2      ka2  0



D0

Estimated values for the model parameters in equations (1)(6), considering the corresponding glucose threshold, are shown in Table 2. Variability of parameters β1 and β2 between subjects is due to a significant difference in peak values of adrenaline in each subject (208.56±116.80 L·min-1 and 2.70±1.81 L·min-1 respectively), although the adrenaline dynamics remains the same among subjects. It thus reflects a change in the system’s gain. Parameters ka1 and ka2 (3.87±1.64 min-1 and 0.07±0.04 min-1 respectively) had an order of magnitude lower than the previous two. Inter-patient variability was lower too reflecting more uniform dynamics among patients.

(15)

(16)

Laplace Transform of the output is: y ( s, )  H1 ( s, )·u ( s )  H 2 ( s, )·x0 ( )

(17)

H1 ( s, )  C ( )[ sI  A( )]1·B ( )  D( )

(18)

1

H 2 ( s, )  C ( )[ sI  A( )] ; x0 ( )  0

As already mentioned, the accuracy of data was tested by the coefficient of determination (R2) and NMRSE. The median coefficient of determination across subjects was 95.45% with a maximum of 99.67% and a minimum 76.22%. The most unfavourable case corresponded to an “outlier” behaviour with a fast increase of glucose level and slow stop of adrenaline secretion which requires further investigation.

(19)

   2·ka 2 1·ka21 ·u ( s ) (20) y ( s, )    2 V ·( k s )( k s ) V ( k s ) ( k s )     a2 e a1 e  

216

9th IFAC BMS Aug. 31 - Sept. 2, 2015. Berlin, Germany Vanessa Moscardó et al. / IFAC-PapersOnLine 48-20 (2015) 213–218

Nine subjects had R2 above 98%, while only three subjects had coefficients of determination below 90%.

217

caused by the oscillatory behaviour of basal secretion, which was neglected. Further investigation is required..

Table 2. Individual estimation of β1, β2, ka1, ka2 and ke, # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

β1

β2

ka1

ka2

ke

(L·min-1)

(L·min-1)

(min-1)

(10-2min-1)

(min-1)

150.93 (0.05) 197.37 (0.01) 119.41 (0.02) 40.35 (0.01) 376.13 (0.01) 197.44 (0.01) 47.06 (0.01) 139.45 (0.01) 172.55 (0.02) 101.41 (0.01) 394.85 (0.02) 161.43 (0.07) 197.46 (0.01) 394.78 (0.01) 177.10 (0.01) 355.26 (0.16) 397.01 (0.01) 294.68 (0.01) 89.05 (0.01) 235.34 (0.02) 160.76 (0.03)

1.46 (0.01) 0.01 (0.04) 2.26 (0.01) 0.01 (0.01) 0.02 (0.01) 2.38 (0.01) 0.01 (0.01) 0.57 (0.01) 5.06 (0.01) 5.02 (0.01) 10.11 (0.01) 0.01 (0.01) 5.06 (0.01) 2.20 (0.01) 0.01 (0.01) 0.01 (0.01) 0.01 (0.01) 12.04 (0.01) 0.40 (0.01) 0.01 (0.00) 10.11 (0.01)

0.52 (0.01) 0.71 (0.01) 5.70 (0.01) 0.09 (0.01) 0.05 (0.00) 0.22 (0.00) 2.33 (0.01) 0.19 (0.00) 1.10 (0.00) 4.66 (0.01) 0.16 (0.00) 1.37 (0.01) 0.03 (0.00) 1.29 (0.00) 4.66 (0.01) 0.56 (0.00) 16.30 (0.0401) 16.30 (0.001) 8.07 (0.001) 16.30 (0.10) 0.75 (0.00)

2.56 (0.01) 8.87 (0.01) 2.33 (0.01) 8.10 (0.01) 1.44 (0.21) 4.25 (0.08) 2.68 (0.01) 2.94 (0.05) 2.05 (0.01) 15.87 (1.02) 12.30 (0.71) 4.56 (0.01) 9.10 (0.11) 18.15 (0.23) 9.30 (0.06) 10.98 (0.01) 6.82 (0.24) 6.45 (0.01) 2.14 (0.01) 3.20 (0.01) 11.67 (0.06)

0.57 (0.00) 0.49 (0.00) 0.12 (0.00) 0.09 (0.00) 2.35 (0.001) 0.81 (0.00) 0.10 (0.00) 0.62 (0.0003) 0.44 (0.01) 0.13 (0.00) 0.53 (0.00) 0.18 (0.0002) 0.20 (0.00) 1.84 (0.01) 0.64 (0.01) 1.25 (0.01) 0.38 (0.00) 0.34 (0.00) 0.02 (0.00) 8.11 (0.01) 6.61 (0.01)

Fig. 3. Model fit in subject #9. (Top) glucose measurements, (middle) input considered in our model, (bottom): adrenaline concentration estimated compared with adrenaline measurements. 3.2 Population parameter estimation The Bayesian analysis provided population estimates of β1pop, β2pop, ka1pop, ka2pop and kepop. Characteristics of them are given in Table 3. Convergence in the Bayesian method was obtained with 12000 samples of the chain. The population mean densities in Figure 4 are narrow and approximately symmetrical. This means that the population mean for the parameters is well defined and its confidence intervals are symmetrical around the mean. Table 3. Population characteristics

The median NRMSE of the model fits was 0.1166±0.067 with a range from 0.045 to 0.284.This indicates a good model fit for all subjects because lower values of it indicate less residual variance. Remark that the different levels of insulinemia (low and high insulin) affected slightly the amplitude values of adrenaline secretion. However, it was not relevant to include an insulindependent factor in the model due mainly to the presence of other more influencing sources of variability.

Population characteristic

Mean(sd)

Interquartile range

β1pop (L-1·min-1)

174.71(28.47)

(153.90-191.40)

β2pop (L-1·min-1)

0.29 (0.24)

(0.14-0.36)

Ka1pop (min-1)

1.15(0.52)

(0.78-1.40)

Ka2pop(min-1)

0.055(0.01)

(0.05-0.06)

Kepop (min-1)

0.49(0.16)

(0.37-0.58)

The population density of β2pop in Figure 4 is strongly skewed to the right, specifying that Type 1 Diabetes patients have low β2 pop, with a most likely value of 0.29 L·min-1, and there is a negligible number of individuals with value of 1 (and more). The same skewed and narrow form is observed in density of ka1 pop. In this case, the likely value is around 1.15 min-1 with no subjects below 0.78 or above 1.40 min-1. On the other hand, values of β1pop are approximately symmetrically distributed around the most likely value of 174.71 L·min-1, with a confidence interval [153.90, 191.40].

In order to better illustrate our model Figure 3 presents a sample (not the best but above average) data fit for subject 9 whose R2 was 97.78% and NRMSE was 0.1127. It is worth remarking that in this example, data showed an increment of adrenaline sooner than the model fit, which may indicate an overestimation of the glucose threshold. This may have been 217

9th IFAC BMS 218 Vanessa Moscardó et al. / IFAC-PapersOnLine 48-20 (2015) 213–218 Aug. 31 - Sept. 2, 2015. Berlin, Germany

It can observed that β1pop and ka2pop density probability is wider than the rest.

Beall, C., Ashford, M. L., & McCrimmon, R. J. (2012). The physiology and pathophysiology of the neural control of the counterregulatory response.American Journal of PhysiologyRegulatory, Integrative and Comparative Physiology, 302(2), R215-R223. Carson, E., & Cobelli, C. (2013). Modeling methodology for physiology and medicine. Newnes. Cryer, P. E., & Gerich, J. E. (1985). Glucose counterregulation, hypoglycemia, and intensive insulin therapy in Diabetes mellitus. The New England journal of medicine, 313(4), 232241. Cryer, P. E. (2001). Hypoglycemia-associated autonomic failure in Diabetes.American Journal of Physiology-Endocrinology And Metabolism, 281(6), E1115-E1121. Cryer, P. (2002). Hypoglycaemia: The limiting factor in the glycaemic management of Type I and Type II Diabetes*. Diabetologia, 45(7), 937-948. Cryer, P. E., Davis, S. N., & Shamoon, H. (2003). Hypoglycemia in Diabetes. Diabetes care, 26(6), 1902-1912. Dalla Man, C., Micheletto, F., Lv, D., Breton, M., Kovatchev, B., & Cobelli, C. (2014). The UVA/PADOVA Type 1 Diabetes Simulator New Features. Journal of Diabetes Science and Technology, 8(1), 26-34. de Galan, B. E., Rietjens, S. J., Tack, C. J., van der Werf, S. P., Sweep, C. G. J., WM Lenders, J., & Smits, P. (2003). Antecedent adrenaline attenuates the responsiveness to but not the release of counterregulatory hormones during subsequent hypoglycemia. The Journal of Clinical Endocrinology & Metabolism, 88(11), 5462-5467. Dejgaard A., Hilsted J. Henriksen JH, Christensen NJ. (1989). Plasma adrenaline kinetics in Type 1 (insulin-dependent) diabetic patients with and without autonomic neuropathy. Diabetologia, 23(11), 810-813. Hansen, N. (2005). The CMA evolution strategy: A tutorial. Vu le, 29. Hansen, N. (2006). The CMA evolution strategy: a comparing review. InTowards a new evolutionary computation (pp. 75102). Springer Berlin Heidelberg. Kovatchev, B. P., Farhy, L. S., Cox, D. J., Straume, M., Yankov, V. I., Gonder-Frederick, L. A., & Clarke, W. L. (1999). Modeling insulin-glucose dynamics during insulin induced hypoglycemia. Evaluation of glucose counterregulation. Computational and Mathematical Methods in Medicine, 1(4), 313-323. Lunn, D. J., Thomas, A., Best, N., & Spiegelhalter, D. (2000). WinBUGS-a Bayesian modelling framework: concepts, structure, and extensibility. Statistics and computing, 10(4), 325-337. Rizza, R. A., Cryer, P. E., & Gerich, J. E. (1979). Role of glucagon, catecholamines, and growth hormone in human glucose counterregulation. Journal of clinical Investigation, 64(1), 6271. Schwartz, N. S., Clutter, W. E., Shah, S. D., & Cryer, P. E. (1987). Glycemic thresholds for activation of glucose counterregulatory systems are higher than the threshold for symptoms. Journal of Clinical Investigation, 79(3), 777-781. Tesfaye, N., & Seaquist, E. R. (2010). Neuroendocrine responses to hypoglycemia. Annals of the New York Academy of Sciences, 1212(1), 12-28. Walter, E. and Pronzato, L. (1997). Identification of parametric models from experimental data. Springer Wilkinson, D. J. (2011). Stochastic modelling for systems biology. CRC press.

When comparing individual results with population characteristic of parameters, it can be realized that density probability of the second is narrower. Interquartile range is reduced too, therefore the variability of parameters is lessened. However, a limitation of this analysis is the reduced number of patients in order to obtain a good representation of the population of patients with Type 1 Diabetes.

Fig. 4. Probability density of parameters β1, β2, ka1, ka2 and ke.

4. CONCLUSIONS This work provides an adrenaline secretion model during hypoglycaemia in patients with Type 1 Diabetes. We consider it as an important step in modelling the complete counterregulation response. The main limitation found has been the variability in adrenaline response, even in the same subject, which may be due to the occurrence of previous hypoglycaemic episodes affecting counterregulatory response. For this, it would be important to validate and corroborate our response by repeating the analysis on different datasets and increasing the number of studied subjects. It would provide subjects with a broader range of adrenaline response variation and a better representation of Type 1 Diabetes population. Future hypoglycaemia detection and prevention systems may improve further with the inclusion of this counterregulatory component.

REFERENCES Amiel, S. A., Sherwin, R. S., Simonson, D. C., & Tamborlane, W. V. (1988). Effect of intensive insulin therapy on glycemic thresholds for counterregulatory hormone release. Diabetes, 37(7), 901-907. Auger, A., & Hansen, N. (2005). A restart CMA evolution strategy with increasing population size. In Evolutionary Computation, 2005. The 2005 IEEE Congress on (Vol. 2, pp. 1769-1776). IEEE

218