فرح محمد ومصطفى Al-Khwarizmi Engineering Journal Al-Khwarizmi Engineering Journal, Vol. 9, No. 1, P.P. 47-59 (2013) Mathematical Modeling of Glucose Regulation System in Term of Perturbed Coefficients Mustaffa Mohammed Basil* Farah Mohammed Ridha Masood** *,**Department of Biomedical Engineering/ Al-Khwarizmi College of Engineering/ University of Baghdad *Email:alfatlaw@msu.edu **Email:fmr82000@yahoo.com (Received 6 February 2012; accepted 19 September 2012) Abstract In this research a recent developed practical modeling technique is applied for the glucose regulation system identification. By using this technique a set of mathematical models is obtained instead of single one to compensate for the loss of information caused by the optimization technique in curve fitting algorithms, the diversity of members inside the single set is interpreted in term of restricted range of its parameters, also a diagnosis criteria is developed for detecting any disorder in the glucose regulation system by investigating the influence of variation of the parameters on the response of the system, this technique is applied in this research practically for 20 cases with association of National Center for Diabetes / Al Mustanseryia University. Keywords: glucose regulation system, compartmental modeling, perturbed coefficients. 1. Introduction In the last decades, there was an increasing demand to obtain a quantitative technique to study physiological systems, for the purposes of diagnosis, therapy, and research. One of the most important techniques that was developed is the compartmental modeling, which is a description to dynamic behavior of physiological systems in term of differential equations based on mass balance equations; these differential equations represent the relationship between exogenous or endogenous material as inputs and the resulted states of physiological system as outputs. The compartmental model was first derived to describe the kinetics of isotopic tracer, science then it was extensively used to deal with wide spectrum of problems in this field, and this can be tracked in [1, 2, 3, and 4]. Compartmental model can be obtained by lumping materials with same characteristics into collections, this will reduce the physiological system into compartments, which can be defined as well mixed and kinetically homogenous materials, and interconnections between them, these interconnections represents the flux of influence from one compartment to another. The recent development in system identification technique has been posed in compartmental modeling. According to system identification theory problem of physiological modeling can be solved through two tasks, first the system specification should be featured by a mathematical model derived from the mass balance equations, secondly the parameters of this mathematical model should be calculated by using the experimental data. The main challenge in the above process is the inherent nonlinearity in physiological system, which is obvious in experimental data; on the other hand a linearized model for these systems exhibits a high degree of uncertainty because of the information loss in the process of linearization. A new algorithm is presented in [5], this algorithm solved this problem by representing nonlinear system by a linear model with perturbed coefficients, and this results a family of models which can cover all aspects of nonlinearity. mailto:Email:alfatlaw@msu.edu mailto:Email:fmr82000@yahoo.com Mustaffa Mohammed Basil Al-Khwarizmi Engineering Journal, Vol. 9, No.1, P.P. 47-59 (2013) 48 Diabetes is one of the most threatening diseases that human face. It is considered one of the major reasons for kidney failure, blindness and limbs amputation. The growth of population and the resulted degradation in health care system, limitation of normal life activity because of the advance in technology, and obesity are the main reasons for the worldwide spreading of this disease, so the importance of formulating a robust mathematical model for the glucose regulation system has been grown over the past decades. This model is derived with the aid of experiments which track the behavior of the glucose concentration after applying an intentioned perturbation in it by specific oral or intravenous dose of glucose. 2. Theory The problem of glucose system regulation in term of perturbed coefficients can be characterized generally as follows: For any experimental set of data obtained from tracking blood glucose concentration after an oral dosage of glucose for fasting person [5]: ( ) i = 1,2, … , N. N: No. of data. let ( )be the nominal function that represents these data, such that ( ) = + +. . + …(1) : number of compartments. , : constant coef icients. there is family of mathematical models: ( ) ∶= ( ): ∈ , , ∈ , , : upper bound of th coef icient. , : lower bound of th coef icient. for = 1,2, … , . The above problem can be solved by an algorithm of three steps: 1. System specification. 2. Nominal parameters estimation. 3. Perturbation ranges calculations. 2.1. System Specification [6] Consider a general physiological system with n compartments as seen in fig (1): Fig. 1. General Compartmental Model. ̇ = − − + …(2) Where: = , = : i exogenous input. : output lux of i compartment. : transferring lux from compartment to . : transferring lux from compartment to . By assuming that this model represents the linearized version of physiological system, then all transferring flux functions can be substituted by linear functions of compartments masses as follows: = …(3) = …(4) = …(5) By substituting eq. (3) to (5) into eq. (2): ̇ = ∑ −∑ − + …(6) After rearranging Eq. (6): ̇ = ∑ − ∑ + + …(7) By assuming that: = − ∑ + …(8) After substituting Eq. (8) in Eq. (7), yields: ̇ = ∑ + …(9) Mustaffa Mohammed Basil Al-Khwarizmi Engineering Journal, Vol. 9, No.1, P.P. 47-59 (2013) 49 The state space representation for n compartments according to the above set of first order differential equations is: ̇ = + …(10) Where: = [ … ] = ⋯ ⋮ ⋱ ⋮ ⋯ = [ … ] Since the experimental data consider the compartments concentration, the nominal output equation is: = …(11) Where: = [ … ] Where: V : Volume of i compartment. The solution for above state space system provides a qualitative understanding for the dynamic specifications of physiological system, so for the following state transition matrix: ( ) = ℒ [ − ] …(12) The time course of compartment concentration as a response to impulse exogenous input can be obtained by: ( ) = [ ( ) (0) + ( ) ] …(13) where: ∶ Amplitude of the impulse exogenous input . eq. (13) results: ( ) = 2.2. Nominal Parameters Estimation [7] For this task, the Least Square Fit algorithm is used, a brief description for this algorithm, that is the process of curve fitting for data set that contains a significant amount of noise and this can be done by minimizing the following function: ( , , … , , ) = ∑ [ ( )− ( )] …(14) Where the weight of experimental data. optimal set of parameters can be obtained by solving the following equation: = 0 = 1,2, … , …(15) The above notation implies that we already have mathematical form of ( ), which was derived previously. 2.3. Perturbation Ranges Calculations [5,8,9] Since least square fit is an optimization technique, then the resulted function doesn't give a full representation to all experimental data, to overcome this weakness in the mathematical model, all unrepresented information will be modeled as a weighted perturbation range for each parameter, so each parameter will be limited by upper and lower bound, this problem can be formulated as follows: ∈ , ∈ , where: = + …(16) = − …(17) = + …(18) = − …(19) , : wights of perturbation. , : lower limits of perturbation. , : upper limits of perturbation. for = 1,2, … , . For mathematical convenience, weights will be represented as: ω = [ω ω ω … … … … ω ] × ω = [ω ω ω … … … … ω ] × While perturbations will be represented as: ̅ = [ … … … … ] × ̅ = [ … … … … ] × Mustaffa Mohammed Basil Al-Khwarizmi Engineering Journal, Vol. 9, No.1, P.P. 47-59 (2013) 50 a) Weights Selection The problem of finding appropriate weight is considered in this part. Weight selection is extremely important for minimizing the family of models, by eliminating unnecessary members. Each parameter in the nominal function has its particular weight, which is defined as the average value of deflections, which occur in nominal function at each time of experimental data, caused by small variation in that parameter. ∆ represents the error between nominal function and experimental data. ∆ = ( )− ( ) …(20) for = 1,2, … , . Sensitivity of nominal function to particular parameter can be defined as: = ( ) …(21) = ( ) …(22) for = 1,2, … , . Now the participation of each parameter in the error between the nominal function and experimental data can be calculated from: ∆ = ∆ ∗ …(23) ∆ = ∆ ∗ …(24) = ∑ + ∑ = ∑ + ∑ for = 1,2, … , . for = 1,2, … , . ∆ , ∆ : the De lection in the nominal function caused by variation in parameter or at each experimental data . For mathematical simplicity, matrix notation will be used in calculation of weights, so first we construct , and ∆ are constructed as follows: = ⋯ ⋮ ⋱ ⋮ ⋯ × = ⋯ ⋮ ⋱ ⋮ ⋯ × ∆ = [ ∆ ∆ … ∆ ∆ ] × The perturbation in nominal function which is de ined as: ∆ = [∆ ∆ ∆ … … ∆ ] × ∆ = [∆ ∆ ∆ … … ∆ ] × and can be calculated as follows: ∆ = ∆ ∗ …(25) ∆ = ∆ ∗ …(26) Finally the weight of perturbation is: ω = ∗ ∆ …(27) = [ … … ] × ω = ∗ ∆ …(28) = [ … … ] × b) Parameter Interval Identification In this part the range of parameter perturbation is calculated, and this can be done by solving the following equation for the variable & at each time of experimental data: ( ) = ( + ) ( ) + ⋯ + ( + ) ( ) …(29) The solution of above n variables can be estimated according to the following: Lets de ine the ℓ range of perturbation in ̅ ∶ ̅ = 0, ≠ ℓ ̂ ℓ, = ℓ …(30) ̅ = 0, ≠ ℓ ̂ℓ, = ℓ … (31) for = 1,2, … , . by substituting eq. (30)in(29), yields ∶ ( ) = + ⋯ + ℓ ℓ + ( ℓ + ℓ ̂ℓ) ℓ + ℓ ℓ + ⋯ + …(32) by substituting eq. (31)in(29), yields ∶ Mustaffa Mohammed Basil Al-Khwarizmi Engineering Journal, Vol. 9, No.1, P.P. 47-59 (2013) 51 ( ) = + ⋯ + ℓ ℓ + ℓ ( ℓ ℓ ℓ) + ℓ ℓ + ⋯ + …(33) Solving eq. (31) and (33) yields ̂ ℓ and ̂ ℓ, which represents the maximum deflection in ℓ parameter caused by the difference between experimental data and nominal function. The actual value of deflection can be calculated by using the following equation: ℓ = ̂ ℓ ∗ ℓ∑ ∑ …(34) ℓ = ̂ ℓ ∗ ℓ∑ ∑ …(35) for = 1,2, … , . Now, for( ℓ = 1,2, … , ), we will have ∶ ̅ℓ = [ … … ] × ̅ℓ = [ … … ] × for upper and lower range of perturbation for parameters respectively = max { 0, } = min { 0, } = max { 0, } = min { 0, } 3. Experimental Results for Oral Glucose Tolerance Test In order to maintain glucose concentration in an adequate level, two hormones are functioning opposite to each other inside the human body, the first hormone is Glucagon, which is produced by a pancreatic islets called Alpha Cells, the secreting of this hormone stimulate the liver to convert the stored glycogen into glucose and release it to the blood, on the other hand the second hormone which is called Insulin, produced by a pancreatic islets called Beta Cells, this hormone assist in assimilation of the glucose. The best approach to obtain data for glucose regulation system modeling, is to observe time response of glucose utilization in blood, one of the ways to do that is the Oral Glucose Tolerance Test (OGTT), it is based on applying oral dosage of glucose, which can be modeled mathematically as impulse function, then the glucose concentration in blood with measuring is measured appropriate sampling rate. Since only the blood glucose concentration will be taken in consideration, then one compartment can represent the system here, and the best function that can be fitted to the experimental data is: ( ) = …(36) In this research, OGTT test is carried out for a 20 person, in association with the National Center for Diabetes / Al Mustanseryia University. The results of the OGTT test base on 75 of oral glucose are exhibited in Table(1). After substituting resulted values (for each case in the table above) in eq. (43), the result will be a family of mathematical models for each case. The family members are bounded by the following two boarder models: ( ) = …(37) ( ) = e …(38) < < Where: ∶ upper limit of the family of models. : lower limit of the family of models. The above argument can be proved practically in Fig. (2) to (21), which represent the time response of glucose regulation system for the 20 cases, in these figures, the dotted blue line represents the experimental data, while the upper red line represents the upper bound in Eq. (37), the lower red line represents the lower bound in Eq. (38), and the middle red line represents the nominal values in Eq. (36). Mustaffa Mohammed Basil Al-Khwarizmi Engineering Journal, Vol. 9, No.1, P.P. 47-59 (2013) 52 Table1, 4. Discussions Since the least square fit is an optimization technique, there is a certain level of information loss which can be useful in describing the dynamic of the system. In the control point of view the process of evaluating any system take into consideration two main characteristics, stability and performance. a. Stability: The glucose regulation system with parametric uncertainty is stable if and only if all the members of this family of models are stable, this means that for stability: < 0 …(39) According to Eq. (36), the stability mainly counts on the values of parameter , any negative value in its range, indicates a possibility of instability in the glucose regulation system. In the results in section 3, there are seven cases with such possibility which are cases 2, 3, 11, 15, 17, 19, and 20. If equal weights are assigned for the possibilities of the values of , the possibility of instability for glucose regulation system % for all cases which exhibit negative values of can be calculated according to the following equation: % = ∗ 100% …(40) Table (2) shows the possibility of instability for all tested cases: No. 1 185.966 185.970 173.702 0.00550 0.0150 0.00260 2 168.084 168.086 164.000 0.00220 0.00430 -0.00210 3 167.437 170.000 167.436 0.00117 0.00214 -0.00022 4 188.748 196.000 188.746 0.00308 0.00658 0.00155 5 187.855 192.338 187.846 0.00517 0.0159 0.00235 6 168.276 173.000 168.2747 0.00480 0.00578 0.00367 7 180.478 180.479 180.000 0.00327 0.00448 0.000820 8 189.348 189.3486 185.978 0.000900 0.00155 0.000159 9 200.819 200.818 184.387 0.00421 0.00800 0.00328 10 164.346 164.3465 161.000 0.00460 0.00512 0.00389 11 150.763 168.000 150.7619 0.000793 0.00730 -0.000687 12 222.484 222.486 208.982 0.00776 0.0106 0.00457 13 152.920 154.5770 152.9182 0.00232 0.00494 0.00138 14 171.227 173.715 171.2226 0.00351 0.00794 0.00185 15 142.161 142.164 142.000 0.00340 0.00655 -0.00112 16 162.767 168.000 162.7656 0.00255 0.00770 0.000780 17 134.7700 137.000 134.7680 0.00187 0.00358 -0.000607 18 188.458 193.161 188.456 0.00425 2.7964 0.00269 19 150.430 155.5115 150.4245 0.00289 0.00922 -0.000391 20 176.8654 176.8686 168.7136 0.00214 0.00638 -0.002122 Mustaffa Mohammed Basil Al-Khwarizmi Engineering Journal, Vol. 9, No.1, P.P. 47-59 (2013) 53 Table 2, b. Performance: fast reaction to any sudden increase in blood glucose level is adopted in this research to evaluate the performance of the glucose regulation system, this means that the performance mainly depends on the decay rate . The normal results for the OGTT adopted by the National Center for Diabetes / Al Mustanseryia University mention in Table (3). By using the recommended values in the above table and solving eq. (36) for the value of ′: ′ = − ′ ln ′ for > Y′ …(41) Where: : initial experimental value. ′: recommended decay rate for normality. So for normality: > ′ …(42) Table (4) shows initial values of experimental glucose concentration, the recommended value for each tested case, the experimental value for the decay rate, and the status of each case according to the condition in eq. (42). Case No. Initial Value ( ) Recommended Decay Rate Experimental Decay Rate Status 1 174 0.001811774 0.0026 normal 2 164 0.001318533 -0.0021 abnormal 3 170 0.001617967 -0.00022 Abnormal 4 196 0.002803935 0.00155 Abnormal 5 196 0.002803935 0.00235 Abnormal 6 173 0.001763743 0.00367 Normal 7 180 0.002094287 0.00082 Abnormal 8 186 0.002367535 0.000159 Abnormal 9 190 0.002544847 0.00328 Normal 10 161 0.001164683 0.00389 Normal 11 168 0.001519346 -0.000687 Abnormal 12 209 0.003339099 0.00457 Normal 13 155 0.000848189 0.00138 Normal 14 174 0.001811774 0.00185 Normal 15 142 0.000118205 -0.00112 Abnormal 16 168 0.001519346 0.00078 Abnormal 17 137 0.001040165 -0.000607 Abnormal 18 186 0.002367535 0.00269 Normal 19 156 0.00090178 -0.000391 abnormal 20 170 0.001617967 -0.002122 abnormal Mustaffa Mohammed Basil Al-Khwarizmi Engineering Journal, Vol. 9, No.1, P.P. 47-59 (2013) 54 Table 3, Table 4, 5. Conclusions In the control point of view the process of evaluating any system take into consideration two main characteristics, stability and performance, according to these two main points a new diagnosis criteria can be concluded from the results for detecting any disorder in glucose regulation system, based on the range of perturbation for parameters in its mathematical model. Fig. 2. The Time Response of Glucose -Case (1). Fig. 3. The Time Response of Glucose -Case (2). Fig. 4. The Time Response of Glucose -Case (3). Case no. % Case no. % 1 0.00% 11 8.60% 2 32.81% 12 0.00% 3 9.32% 13 0.00% 4 0.00% 14 0.00% 5 0.00% 15 14.60% 6 0.00% 16 0.00% 7 0.00% 17 14.50% 8 0.00% 18 0.00% 9 0.00% 19 4.07% 10 0.00% 20 24.96% ′( ) ′( ) 60 <200 120 <140 Mustaffa Mohammed Basil Al-Khwarizmi Engineering Journal, Vol. 9, No.1, P.P. 47-59 (2013) 55 Fig. 5. The Time Response of Glucose -Case (4). Fig. 6.The Time Response of Glucose -Case (5). Fig. 7. The Time Response of Glucose -Case (6). Fig. 8. The Time Response of Glucose -Case (7). Fig. 9. The Time Response of Glucose -Case (8). Fig. 10. The Time Response of Glucose -Case (9). Mustaffa Mohammed Basil Al-Khwarizmi Engineering Journal, Vol. 9, No.1, P.P. 47-59 (2013) 56 Fig. 11. The Time Response of Glucose -Case (10). Fig. 12. The Time Response of Glucose -Case (11). Fig. 13. The Time Response of Glucose -Case (12). Fig. 14. The Time Response of Glucose -Case (13). Fig. 15. The Time Response of Glucose -Case (14). Fig. 16. The Time Response of Glucose -Case (15). Mustaffa Mohammed Basil Al-Khwarizmi Engineering Journal, Vol. 9, No.1, P.P. 47-59 (2013) 57 Fig. 17. The Time Response of Glucose -Case (16). Fig. 18. The Time Response of Glucose -Case (17). Fig. 19. The Time Response of Glucose -Case (18). Fig. 20. The Time Response of Glucose -Case (19). Fig. 21. The Time Response of Glucose -Case (20). Mustaffa Mohammed Basil Al-Khwarizmi Engineering Journal, Vol. 9, No.1, P.P. 47-59 (2013) 58 6. References [1] Carson,E.,R.,Cobelli,C.,and Finkelstien,L. "The Mathematical Modelling Of Metabolic And Endocrine System". John Wiley Sons, New York. ( 1983). [2] Godfrey, K.,. "Compartmental Models and their Application". AcademicPress, London,( 1983). [3] Jacquez, J.A.. "Compartmental Analysis in Biology and Medicine". 3rd ed., Biomedware, Ann Arbor, MI,( 1996). [4] Cobelli, C., Foster, D.M., and Toffolo, G. "Tracer Kinetics in Biomedical Research : From Data to Model". Plenum Publishing Corp.,(2000). [5] Basil, M., M., "System Identification Algorithm for Systems with interval Coefficients",under press. [6] Bronzino, J., D., "The Biomedial Engineering Handbook", IEEE Press, (2006). [7] Kivsalaas, J., "Numerical Method in Engineering", Cambridge University Press, (2010). [8] Sanchez, R., S., "Robust System Theory and Application", John Wiley & sons (1998). [9] Bhattacharyya, S., P., Chapellat, H., and Keel, L., H., "Robust Control Parametric A p p r o a c h " , P r e n t i c e H a l l ( 1 9 9 5 ) . )2013(47-59 ، صفحة1، العدد9مجلة الخوارزمي الھندسیة المجلدمصطفى محمد باسل 59 النمذجة الریاضیة لمنظومة السیطرة على سكر الكلوكوز بصیغة المعامالت المتذبذبة **فرح محمد رضا مسعود* مصطفى محمد باسل جامعة بغداد /كلیة الھندسة الخوارزمي / قسم ھندسة الطب الحیاتي** ،* alfatlaw@msu.edu :البرید االلكتروني* fmr82000@yahoo.com :البرید االلكتروني** الخالصة باستخدام ھذه الطریقة ،في ھذا البحث تم استخدام طریقة مستحدثة للنمذجة بصورة عملیة في عملیة نمذجة منظومة السیطرة على تركیز السكر في الدم ریاضي سیتم الحصول على مجموعة من النماذج الریاضیة بدال من نموذج واحد لتعویض الخسارة فى المعلومات الناجمة عن عملیة استحصال اقرب نموذج ت في النموذج الریاضي بین قیمتین ھذا التنوع في النماذج الریاضیة داخل المنظومة الواحدة ناجم عن تغییر المعامال. یحاكي المنظومة بصورة مقربة ایضا في ھذا البحث تم تطویر طریقة تشخیصیة ریاضیة الكتشاف اي خلل في منظومة السیطرة على السكر في الدم من خالل مراقبة تاثیر ھذا ،محددتین لة بالتعاون مع المركز الوطني للسكري فى الجامعة و تم تطبیق ھذه الطریقة بصورة عملیة على عشرین حا. التغییر في المعامالت على استجابة المنظومة .المستنصریة mailto:alfatlaw@msu.edu mailto:fmr82000@yahoo.com