Table Of Content

Journal of Information, Control and Management Systems, Vol. 4, (2006), No.2 203 MODEL OF A PLASMA RENIN ACTIVITY AFTER NICARDIPINE TREATMENT Anna TOLEKOVA, Kaloyan YANKOV Thracian University, Medical Faculty Stara Zagora, Bulgaria e-mail: [email protected], [email protected] Abstract Mathematical model of plasma renin activity after nicardipine treatment is developed. System identification of the process is done applying the cyclic coordinate descent as optimization procedure. The model allows to predict the effects of different drug doses and allows the researcher to examine the behavior of the system under all conceivable conditions. Keywords: system identification, mathematical model, cyclic coordinate descent, renin angeotensin system, plasma renin activity. 1 INTRODUCTION The values of blood pressure are an important indicator of human cardio-vascular status. A considerable part of the long-term regulation is exerted by the renin- angiotensin system. The deviations of the blood pressure out of reference limits are treated with different groups of drugs, that influence various regulatory factors, including plasma renin activity (PRA). The later is an integral index of renin- angiotensin system state. The changes of PRA and their regulatory kinetics have been studied after physiological stimulation [1] and after treatment with different drugs: inhibitors of calmodulin [2], blockers of prostaglandin syntase [3], angiotensin-converting enzyme inhibitors [4] and calcium channel blockers [5, 6] on male Wistar rats. Their dynamic characteristics are investigated employing statistics methods and numerical analysis [7, 8, 9, 10]. The received data are insufficient and it is incorrect to translate the results and the conclusions upon humans. In general resources and methods are necessary in order to predict the effects of certain treatment on humans and then to pass to the clinical research phase of that same treatment. An alternative approach, which proves its effectiveness for system evaluation and determination the quality of regulation at reasonable price, is the use of mathematical models [11]. The power of modeling lies in its abstraction, since a single family of models may describe many real systems. Mathematical modeling means to find out the relations that characterize the internal 204 Model of a Plasma Renin Activity structure of the system and to describe the interdependencies between the input and output variables when the system is aftected by the environment in a particular way. The purpose of this work is a mathematical model formulation of PRA dynamics under the influence of different doses of nicardipine. The model will be used for investigation of PRA from the point of view of system theory. 2 EXPERIMENTAL DESIGN The experiments were carried out on 145 male white Wistar rats. PRA was assessed radioimmunologicaly (DiaSorin-Biomedica Ltd.) after perorally application of Nicardipine in doses 10, 20, 40 60 mg/kg body weight (b.w.). PRA was sampled in the intervals showed in a Table 1. Animals received humane care and the study complied with the Institution's guidelines of Thracian University and with the Guidelines for Breeding and Care of Laboratory Animals (Veterinary Public Health Reports, 1994). Time Dose [mg/kg] b.w. [hours] 10 20 40 60 0 7.58±0.8 7.58±0.8 7.58±0.8 7.58±0.8 0.5 27.13±3.1 31.68±3.1 32.13±2.9 33.80±6.8 1 35.18±4.5 40.92±9.4 44.44±1.7 51.22±11.2 3 40.07±6.4 45.91±3.4 48.56±2.6 49.61±8.7 5 23.21±4.3 27.55±4.5 28.84±3.9 31.18±5.4 7 11.17±2.8 12.47±1.5 14.58±1.5 17.24±3.8 9 7.58±0.8 7.83±0.9 8.85±0.9 9.80±2.1 11 7.56±0.6 7.50±0.7 7.80±0.6 7.58±0.4 Table 1. Plasma renin activity [ng/ml/h], presented as means and standard deviation. 3 DESIGN OF SYSTEM IDENTIFICATION PROCEDURE An identification experiment is performed by exciting the system using some sort of input signals and observing a change in one of the system state variables as an output signal over a time interval. For modeling of the PRA production process, the experiment was designed to observe the following sequence: • Input signal U(t). A short perorally application of nicardipine is considered as Dirac function. The signal amplitude is correlated to the drug dose. • Output response y(t). During the experiment a few values of the output signal PRA are recorded. Let Φ(t) ⊂ y(t) be the observational vector, which is known from the experiment: Φ(t) = [ϕ , ϕ ,... ϕ ,]Т, where N is the number of samples 1 2 N The vector Φ(t) is used during the identification process. Journal of Information, Control and Management Systems, Vol. 4, (2006), No.2 205 • Identification time t . That is the time for process modeling. After this time the p system response reaches the steady state level and then the system state variables are time independent. The maximal duration time was fixed to 11 hours. • Sampling time. The first two samples were taken at 30 min and 1 hr and the subsequent were taken at every 2 hrs (see Table 1). • Structure of the proces model: Consists of the mathematical equations which describe the behaviour of the process. The precision of the model depends on the appropriate choice of equations. Let y(t,q ,q ,…,q ) is the model for identification, 1 2 m t is the time, Q = [q ,q ,…,q ] is the identification vector make of the 1 2 m parametersin terms of which we want to describe the system behavior. Defining the values of those parameters is the object of the system identification. PRA follows an oscillation curve. The most appropriate model is the second order ordinary differential equation: d2y(t) dy(t) +2ζϖ +ϖ2y(t)+K = KϖU(t) (1) dt2 dt 0 u where: U(t) – the input signal. In our case the applied dose of nicardipine. ζ(d) – the damping ratio. ω(d) - the undamped, natural frequency of the system. К (d) – the base level. 0 K (d) – the sensibility of the process to the input influence (proportionality u coefficient). The parameters above are unknown and they must be calculated in order to identify the process. All of them are dose (d) dependent and they form the identification vector Q(d): Q(d)=Q(ζ(d), ω(d), К (d), K (d)) 0 u • Identification method. The mathematical model is determined using some identification methods. For nonlinear models very few results have been obtained and no standard algorithm exists for testing global identifiability. Various approaches have been proposed, e.g., power series [12], differential algebra [13], stochastic approximation [14], similarity transformation methods [15]. As an identification method we use cyclic coordinate descent (CCD) [16] realized by software KORELIA-DYNAMICS [17]. The values ϕ of the model y(t,Q) are i known in m number of points t : i y(t, q , q ,…,q ) = ϕ(t), i=1,2,…,m (2) i 1 2 m i i The identification vector Q must be defined from the equation system (2). But this system is not well defined. For each point ϕ we define the residual d: i i d(Q) = ||y-ϕ|| i i i The aim is to minimize d for the identification time t : i p 206 Model of a Plasma Renin Activity D(Q) = inf|| d(Q) || < ε, ε>0 (3) i Thus the identification goal is translated into an optimization problem. The CCD method is an iterative heuristic search technique that attempts to minimize D(Q) by varying one parameter at a time. It reduces the multi-criteria optimization to the single- criteria one. A number of iterations are made over the model equation to find the global minimum D for the system (3). The successive approximation of the identification parameters q on the k-th iteration is evaluated using the procedure: i qk = argminD(qk,qk,...θ,qk−1,...,qk−1) i 1 2 i i+1 N (4) θ where argmin(…θ…) is the value of θ for which D(Q) has a minimum on the corresponding coordinate while the other coordinates are fixed. The implemented criterion to improve the estimation of the model parameters is the uniform fitting: D(Q)=max|d |→min (5) i The criteria for termination of the identification process is the unequality: | D(Qк) - D(Qк-1) | < δ, where δ is an accuracy estimation (6) CCD is simple for realization, relatively robust method. The main disadvantage is its locality and the fact that it may take a lot of calculations. 4 REZULTS The system identification is made using as a model the differential equation (1) at initial conditions: dy(0) y(0) = 7.58 =0 (7) dt The calculated values of the ζ(d), ω(d), К (d) and K (d) using CCD are presented in 0 u Table 2. Dose [mg/kg] b.w. parameter 10 20 40 60 ζ(d) 0.67 0.95 1.35 1.58 ω(d) 0.64 0. 64 0.64 0.64 К (d) -5.97 -3.00 -1.62 0.42 0 K (d) 30.00 25.83 19.03 15.18 u Table 2. Calculated values of the identification parameters 4.1 Damping ratio ζ(d). The dose-dependent values are graphically showed on Figure 1. The graph could be described with the exponential growth model equation of the type : d+∆ − (8) ζ(d)=C (1−e D )+C ∞ const Journal of Information, Control and Management Systems, Vol. 4, (2006), No.2 207 Where the unknown parameters for identification are: C = ζ(d→∞) ∞ D – dose-constant. ∆ – dose correction parameter; C – free term const Using CCD, the calculated values are: C =6.6; D=123.4; ∆=73.9; C = -2.6 ∞ const and the dose-dependent equation for the damping ratio is: d+73.9 d+73.9 − − (9) ζ(d)=6.6(1−e 123.4 )−2.6=4−e 123.4 Figure 1. Damping ratio ζ(d) 4.2 Proportionality coefficient K (d). u Figure 2. Proportionality coefficient K (d) u 208 Model of a Plasma Renin Activity This is an exponential decrease model of type: d+∆ − (10) K (d)=C e D +C U 0 const Where the unknown parameters are: C = K (0) 0 u D – dose-constant. ∆ - dose correction parameter; C – free term const Using CCD the calculated values are: C =24.28; D=47.6; ∆=-457, and the ∞ obtained equation is: d−457 − (11) K (d)=24.28e 47.6 u 4.3 Constant base level К (d) – figure D. 0 Some models were approved and the equation that gives minimal residual is: K = C *ln(d ) + C (12) 0 0 const After optimization the calculated coefficients leads to the equation: K (d)=3.64ln(d)−14.47 (13) 0 Figure 3. Constant base level K 0 Finally, the time and dose dependent plasma renin activity model is described by the system of equations: d+73.9 − ζ(d)=4−e 123.4 −d−457 (14) K (d)=24.28e 47.6 u K (d)=3.64ln(d)−14.47 0 Journal of Information, Control and Management Systems, Vol. 4, (2006), No.2 209 d2y(t,d) dy(t) +1.27ζ(d) +0.406y(t)+K (d)= K (d)U(t) dt2 dt 0 u With initial conditions: dy(0) y(0) = 7.58 =0 dt The graphics of the experimental data interpolated using cubic spline and generated models of PRA for dose of 10 mg.kg and 60 mg/kd are shown on Figure 4. Figure 4. Experimental data and simulation curves of PRA after treatement with doses 10 and 60 mg/kg b.w. 5 TRANSFER FUNCTION The transfer function summarizes a system in the Laplace domain. It is the ratio of the output signal to the input one. A general second-order transfer function is: K ω2 G(s)= U (15) s2 +2ζωs+ω2 After substituting the corresponding coefficients from system (14) we obtain the dose-dependent transfer function: d−457 − 9.9451e 47.6 G(s,d)= d +73.9 (16) − 2 123.4 s +1.28(1 -6.6e )s+0.4096