Total Mass TCI driven by Parametric Estimation

Total Mass TCI driven by Parametric Estimation Margarida M. Silva∗, Cl´ audia Sousa†, Raquel Sebasti˜ ao‡, Jo˜ ao Gama§, Teresa Mendonc¸a¶, Paula Rochakand Sim˜ ao Esteves∗∗ ∗DMA - Faculdade de Ciˆ encias da Universidade do Porto (FCUP), Porto, Portugal, Email: margarida.silv[email protected] †Escola Superior de Educac¸˜ ao Jean Piaget - Instituto Piaget, Vila Nova de Gaia, Portugal, Email:[email protected]g ‡LIAAD - INESC Porto, L.A. and FCUP, Porto, Portugal, Email: [email protected] §LIAAD - INESC Porto, L.A. and Faculdade de Economia da Universidade do Porto, Porto, Portugal, Email: [email protected] ¶DMA - FCUP, Porto, Portugal, Email: [email protected] kFaculdade de Engenharia da Universidade do Porto, Porto, Portugal, Email: [email protected] ∗∗Anesthesiology Department - Hospital Geral de Santo Ant´ onio, Porto, Portugal, Email: simao.este[email protected] Abstract—This paper presents an algorithm, Total Mass Target Controlled Infusion (TCI) system, including n on line parameter estimation scheem along with a a mass control for compartmental systems. The system comprises an On line tuned Algorithm for Recovery Detection (OLARD) applied after an initial bolus administration and a Bayesian identification method based on sparse measurements on the accessible signal which aims to minimize a Bayesian cost function. To design the drug dosage profile, two algorithms are here proposed. During the transient phase, an Input Variance Control (IVC) algorithm is used. It is based on the concept of TCI and aims to take the drug effect to a predefined target within an apriori defined interval of time. Afterward the drug dose regimen is controlled by a Total Mass Control (TMC) algorithm. The mass control for compartmental systems is robust even in the presence of parameter uncertainties . The whole system feasibility have been evaluated for the case of Neuromuscular Blockade (NMB) level and was tested both in simulation and in real cases. I. INTRODUCTION Automation in biomedicine undoubtly plays an important role in the improvement of health care, which claims for robust and secure technologies. This is particularly true when we want to determine and apply the adequate drug dosage regimen to each patient. In general anaesthesia environment, the major difficulty concerning the drug dose administration is the high level of uncertainty in the model relating the manipulated variable (amount of drug) and the choice of the adequate sensor our even process variable (for instance the neuromuscular blockade level, NMB, or the depth of anaesthesia evaluated by brain bispectral index of unconsciousness, BIS). Although the structure of the pharmacokinetic/pharmacodynamic (PK/PD) models is relatively well known, there is a high parameter variability from patient to patient that suggests the use of adaptive methods both for identification and control purposes. Trying to solve the patient identification problem and the control of drug dose administration, several approaches have recently being developed [1], [2]. Actually, Target Controlled Infusion (TCI) strategies are commonly accepted and used due to its specific features. TCI allows for a controlled infusion in such a manner as to attempt to achieve a user-defined drug concentration in a body compartment or tissue of interest [3]. In order to reach that PD concentration, the dose profile is predicted by validated population models. Therefore, these open-loop control devices do not compensate for a mismatch between the models and the patient dynamics and they are not completely adequate for drug administration [4]. This study presents the development of an integrated control algorithm that deal with the need for adaptation in TCI. The proposed control strategy computes the adequate drug dosage regimen that drives the drug effect to a desired target in a prespecified period of time and determines the loading dose that should be administered to control the system around that working point. The global algorithm architecture also incorporates a identification method based on a Bayesian optimization for parameter identification and a online tuned algorithm for recovery detection [5]. The implementation of the whole system controller is here illustrated with an example of closed-loop control of Neuromuscular Blockade (NMB) with the muscle relaxant atracurium. The remainder of the paper is structured as follows: the NMB model structure together with the identification algorithm are defined in Sections II-A and II-B. The online tuned algorithm for recovery detection is summarized in Section II-C. The input variance control algorithm and total-mass control algorithm are explained in detail in Sections II-D and II-E. Some simulation results for the NMB case-study are given and discussed in Section III, also with one real case where the system has been tested. A brief summary of the main proposed advances concludes the paper (Section IV). II. ADAPTIVE TCI STRATEGY A. Neuro Muscular Blockade Model The dynamic response of NMB for several muscle relaxants may be modeled by a Wiener structure: a linear part followed by a non-linear static dynamic one. For the muscle relaxant atracurium the resulting equations are the following. The linear PK part relates the drug infusion rate u(t)[µg kg−1min−1] with the plasma concentration cp(t) [µg ml−1] (1), cp(t) = u0X i=1,2 aie−λit(1) A single bolus of u0, given at t= 0 min (in the induction phase) may be described by u(t) = u0δ(t)µg kg−1. For atracurium, the typical bolus is of 500 µg kg−1. The PD has two linear parts relating cp(t)and the effect concentration ce(t)[µg ml−1](2), and a non-linear static Hill relationship between ce(t)and the blockade level r(t)[%] (3). The vactor Θincorporates the patient dependent parameters: Θ = {a1[kg ml−1], λ1[min−1], a2[kg ml−1], λ2[min−1], λ[min−1], τ [min−1], C50 [µg ml−1], γ (dimensionless)}. ce(t) = λ 1−λτ e−λt −λ 1−λτ e−(1/τ)tcp(t)(2) r(t) = 100Cγ 50 Cγ 50 +cγ e(t)(3) Remark that in practice the intermediate signal ce(t)is not measured and its theoretical value is used. The NMB level of a patient is typically quantified between 0% (full paralysis) and 100% (full muscular activity). For control purposes, during the period where the bolus is acting, the value of the reference is fixed at a low level, being gradually raised to the set-point r∗(tipically 10%) in order to avoid sudden changes, stabilizing on this value after minute 75 (steady-state) [6]. Taking into account this parameterization and in order to cover a wide range of behaviours, a bank of nonlinear dynamic models M={Mi(Θi P K ,Θi P D)i=1,...,100}was generated using the probabilistic model discussed in [7]. Since the average parameters do not convey information regarding the spread of the individual values, it is reasonable that the use of additional information on the interindividual variability of the PK/PD parameters may lead to techniques achieving improved results. B. PK/PD Bayesian Identification When a measurement of the response signal r(t)is available, the estimates of the patient PK/PD parameters are actualized by minimizing the following Bayesian objective cost function, which assumes an underlying log-normal distribution for Θ[8], H(ˆ θ, u) = m X j=1 (ln(θj)−ln(ˆ θj))2 σ2 j +(4) + n X k=1 (ln(T(tk)) −ln(T(tk;u, ˆ Θ))2 σ2 k with T(tk) = ln r(tk) 100 −r(tk)=γln C50 ce(tk)(5) where the vectors Θ=[θ1, θ2,· · · , θm]and ˆ Θ=[ˆ θ1,ˆ θ2,· · · ,ˆ θm], (with m= 7,since for C50 was assumed the mean population value (??)), represent the PK/PD for some population model and for the individual patient, respectively; u(t)is the dosage regimen; T(tk)is a function of the NMB measured signal at time tk;T(tk;u, ˆ Θ) is the predicted at time tkinduced by u(t)on the patient model with parameter vector ˆ Θ;σ2 jis the variance of the logarithmic distribution assumed for the population PK/PD parameter θj; and σ2 kdescribes the variance of the error on the measurement of the individual NMB signal at t=tk. C. The On line tuned Algorithm for Recovery Detection - OLARD Considering an automated drug dose administration, it should be noted that the time evolution of the measured response r(t)is highly dependent on the time instant t0that is chosen to start the continuous drug infusion u(t). Therefore, an on-line detection for the assumed initial recovery t0is of outstanding importance and it may be used as a decision flag, alarm or advisory component in an automatic drug delivery setup. During the induction phase, the proposed On line tuned Algorithm for Recovery Detection (OLARD) initially computes a detection baseline in order to overcome a variety of sensor drawbacks. Thereafter, the detection of the NMB recovery from the initial bolus administration is computed through an empirical algorithm developed taking into account the baseline and the characteristics of the signal under study. The algorithm parameters had been tuned using both in simulation, using the model bank Mand in a collection of real cases previously obtained during monitor or control procedures in patients under surgery. Moreover, in order to achieve a robust and reliable performance, the default algorithm parameters can be adjusted on line, taking into account the predicted degree of patient variability and the level of sensor noise, as well as the clinical experience about the system environment. D. Input Variance Control algorithm - IVC The Input Variance Control algorithm (IVC) was designed to control the system in the transient phase [t0, t1], where t0 is the previously identified time for recovery and t1is the instant time where steady-state is supposed to be reached, i.e., r(t1) = r∗. The IVC algorithm approach designs an optimal dose regimen profile for taking the NMB level from its value r(t0)at t0to a prespecified value r(t1) = r∗at t1. Assuming linear time-invariant dynamics, the PK/PD system drug response may be modeled by the following discretetime state-space system: x(t+ 1) = Φ x(t) + Γ u(t) y(t+ 1) = C x(t+ 1) (6) where u(t)∈Ris the input (piecewise constant drug dose infusion), x(t)∈Rn×1the state-vector, y(t)∈Ris the output (effect concentration of the drug), C∈R1×n,Φ∈Rn×nand Γ∈Rn×1(with nequal to the system order). It should be stressed that Φand Γare dependent on the values of Θ. This dependence on the parameters is omitted only for the sake of simplicity in representation. Given an initial time t0, an initial state x(t0) = X0, and an input signal u(t)defined for all t, a solution of (6) for t > t0 is given by [10]: y(t) = Cφ(t, t0)X0+ t−1 X τ=t0 Cφ(k, τ + 1)Γu(τ)(7) where φ(i, j) = Φi−jfor i≥j+ 1. The evolution of the system from t0(the beginning of drug infusion) to t1=t0+p(the time when the steady-state is supposed to be reached) can then be represented by Y= ΣX0+SU(8) where YT=y(t0+ 1) y(t0+ 2) ··· y(t0+p−1) y(t1) ΣT=CΦCΦ2··· CΦp−1CΦp UT=u(t0)u(t0+ 1) ··· u(t0+p−2) u(t1−1) S=         CΓ 0 0 ··· 0 CΦΓ CΓ 0 ··· 0 . . .. . .. . .. . .0 ··· ··· ··· CΓ 0 CΦp−1ΓCΦp−2Γ··· CΦΓ CΓ         The reference profile in the effect concentration is obtained from the pre-fixed c∗ e(t1) = f−1 NL(Θ, r∗)and c∗ e(t0) = f−1 NL(Θ, r(t0)): c∗T e=c∗ e(t0)c∗ e(t0+ 1) . . . c∗ e(t1)(9) Making Y=c∗ ein (8), and in order to avoid too large control signals and syringe pump saturation, a variance constrained input strategy was developed based on [11]. Having already calculated one initial guess of the steadystate drug dose as being a non linear function on the system parameters and target value in NMB (10), uss =gNL(Θ, r∗) = C50 ((100/r∗)−1)1/γ (a1/λ1) + (a2/λ2)(10) instead of aiming to determine the total amount of drug to give to the patient, we may want to determine which amount of drug must have to be added or decreased to that value for the effect concentration ce(t)to follow the imposed reference profile c∗ e. The PK/PD system description can then be seen as incremental, and u(t)may be decomposed in U=Uss +∆U. Hence, (8) is modified into Y= ΣX0+SUss +S∆U(11) where UT ss =uss uss ... ... uss ∆UT=∆u(t0) ∆u(t0+ 1) ... ∆u(t0+p−1) The adaptive control problem can then be formulated as follows. Problem 1. Given c≥0, find an ∆Usolving the Lagrangean function L(f, ρ) = kS∆U−Dk2+ρ{k∆Uk2−c2}(12) where ρis as a Lagrange multiplier satisfying ρ(k∆Uk2−c2) = 0 and D=Y−ΣX0−SUss. Note that, in each time instant the patient response r(t) is measured by the adequate sensor after applying u(t0) = uss + ∆u(t0)corresponding to the first value of the obtained optimal drug regimen (U=Uss + ∆U). Thereafter, inverting the Hill equation (3), the corresponding effect concentration is calculated and a new target profile c∗ e=Yis considered. This procedure is then repeated, until the end of the transient period is reached, assuming in each iteration la new transient interval [t0+l, t1], l = 1,· · · ,(t1−t0+ 1). E. Total-mass control algorithm - TMC For practical purposes and without loss of generality, the identification procedure is stopped as soon as the parameters are assumed to be reasonably identified and the reference signal attains the predefined value r∗= 10%. It turns out that in order to avoid computational bender, from this time instant the individual patient parameterization is assumed to be obtained. Neverthless, it is possible to on-line require suplementary identification procedure in order to obtain a more accurate parameter estimation. From this point on, the control strategy is changed to a mass control scheme which has proved to have good robustness properties under the presence of uncertainties [12]. For this purpose the patient is modelled as a three-compartmental system. The drug infusion dose u(t)is administered in the central compartment, according to the law, u(m) = max (0,˜u(m)) ˜u(m) = P3 i=1 qimi+λ(M∗−M(m)) (13) where M(m) = Pi=1,2,3miand m= (m1, m2, m3)T.mi is the amount of drug in compartment iand qi≥0is a rate constant such that qimiis the flow rate from compartment ito the environment. The values of qias well as the initial values of mi,i= 1,2,3, are obtained from the previously estimated parameters. M∗is a suitable value for the total amount of drug in the three compartments, computed in such a way that the corresponding effect concentration yields the target value of r∗for neuromuscular blockade [12]. Since the parameters qi’s are affected by an estimation error, a steady state error will occur. Indeed, denoting the error of qiby ∆qi, the previous control law becomes u(m) = max (0,˜u(m)) ˜u(m) = P3 i=1 (qi+ ∆qi)mi+λ(M∗−M(m)) (14) and it is shown in [12] that, when the drug infusion dose u(t) is administered in the central compartment according to this law, the amount of drug in each compartment will converge to a value that depends on ∆qi’s. According to the robustness analysis that was made in [12], the decreasing of the values ∆qi, i = 1,2,3, would correspond to the decreasing of the steady state error. However, ∆q1,∆q2and ∆q3are fixed, after the beginning of this control scheme. Therefore, in order to reduce the steady state error, we redesign the mass control law, by replacing the original value of M∗by another one, computed based on information about the obtained steady state error. A more complete description of this strategy will be given elsewhere. Fig. 1. NMB response (upper plot) and drug dose profile (lower plot) for patient number 60 ( ˆ Θ = Θ60) in Mthat result from the control with the Total Mass TCI. The ’dot’ in the xx-axis indicates time instant t0computed by OLARD. Note that NMB signal is the filtered NMB simulated superimposing noise (σ= 0.3; 0). Fig. 2. NMB response (upper plot) and drug dose profile (lower plot) for patient number 60 ( ˆ Θidentified with the Bayesian algorithm from t0until t1) in Mthat result from the control with the Total Target Control Algorithm. The ’dot’ in the xx-axis indicates time instant t0computed by OLARD. Note that NMB signal is the filtered NMB simulated superimposing noise (σ= 0.3; 0). III. RESULTS For the global approach developed, comprising the dedicated identification and control algorithms previously described, simulation studies have been carried out using the model bank M. Notice that to achieve a high level of NMB (neuromuscular blockade) in a short period of time, in order to fulfill clinical requirements, a typical 500 µg kg−1bolus of atracurium is administered in the beginning of surgery (to the bank of simulated models and real patients). After the administration of the bolus, the label of NMB increases very quickly (the variable rthat measures muscular activity decreases), and full paralysis is induced in a few minutes. Following the initial period, the control objective is to follow a specific reference profile with a target value. Fig. 1 and 2 illustrate the filtered response, the target value and the administered drug dose profile for patient M60 ∈ M after application of the developed Total Mass TCI strategy. In Fig. 1 and in order to highlight the algorithm performance, the exact parameterization of the simulated patient ( ˆ Θ=ΘM60 ) was assumed to be known. The NMB response is represented just to a maximum value of 30% for the sake of illustrate the reference tracking in detail. The dot indicates the beginning of recovery estimated according to OLARD specifications, and the instant time t0is thereafter used as the decision flag to enable the beginning of the continuous infusion of atracurium administration. After the t0identification, the value of NMB is driven, by IVC strategy, from r(t0)in order to achieve the target value r(t1= 75min) = r∗= 10%. Afterwards, the TCM is applied to control the system by designing the adequate drug dose profile to attain the total mass M∗, corresponding to r∗(t) = 10%. Following the TMC strategy the NMB level is driven to a value bellow 10%, exhibiting a negleted constant steady state error, since it was assumed the existence of parameters uncertainties (∆qi6= 0). Fig. 3. Evolution of parameter values, from minute 10 until t1= 75 min, corresponding to the simulation represented in Fig. 2. Refinements of the parameterization were carried out from t0= 30 min until t1. The values are normalized by the respective true parameter values of patient 60 in M. For the case illustrated in Fig. 2, during the transient period [t0, t1], the patient parameters were identified every Fig. 4. Total mass evolution for the simulation control essay with model 60 (ˆ Θ). In solid line is represented the target level in total mass (M∗(ˆ Θfinal)); in dotted line is represented the theoretical value of M∗(Θ60)for this simulated patient; in dashed line is represented the time evolution of the system total mass (M(m)); in dots is represented the different values of M∗(ˆ Θk)in the transient phase 3 minutes by the identification Bayesian algorithm. Time parameter values evolution is shown in Fig. 3 being the values normalized by the true value of each one ( ˆ Θ/ˆ Θ60). Notice that the last two parameter refinement procedures do not introduce relevant additional information about the system dynamics, indicating that the samples obtained from the system do not provide new information about it. Moreover, apart from parameters λ1and λ2that have been accurately identified, the other parameters present an identification error less that 10%. Remark that, in this case, during the identification procedure the majority of the parameters present a smooth trajectory from the initial guess to their final estimates, whereas the parameter a1presents a significant variation (following the typical scenario observed in the model bank M). The NMB level r(t)presents an overshoot after the beginning of continuous infusion, since the drug profile initially calculated through IVC algorithm relies on initial ’rough’ parameter estimates. However, it can be clearly observed that, with the improvement of the parameter estimates during the transient period, the drug dose rate profile is being updated accordingly, driving the NMB value to the desired target. After, TMC strategy sustains the continuous infusion and brings the NMB level to a value near 8%. This negleted constant steady state error, is due to the use of parameter estimates, as referred before. Fig. 4 represents M∗(ˆ Θfinal), the desired value for the total mass of the system (solid line) and the time evolution of (M(m)), the current total mass (in dashed line). At the time instant indicated by the arrow the control law was changed from (13) to the referred redesigned mass control law, by replacing the original value of M∗by another one, computed based on information about the obtained steady state error. The total mass convergence is then modified and the total mass then goes into the target value of M∗. The Total Mass TCI strategy was recently integrated in the software Hippocrates [[13]] and is actually under evaluation on patients submitted to general anaesthesia. Fig. 5 illustrates the first case so far obtained on a patient undergoing elective surgery. The clinical performance was considered satisfactory in spite of the steady-state error observed during the application of TMC phase that was higher than the observed during simulation. The filtered NMB signal indicates the presence of sensor noise superimposed on the measured signal and a quite high baseline. The recovery estimated according to OLARD specifications and denoted by the ’dot’ was considered correct. Time parameter values evolution is shown in Fig. 6 being the values normalized by the initial guess corresponding to the population parametrization true value of each one ( ˆ Θ/¯ Θ). Notice that, as in the reported simulated cases, the last two parameter refinement procedures do not introduce relevant additional information about the system dynamics. In the real case some of the parameters present significant variations from the population parameters. In Fig. 7 is represented M∗(ˆ Θfinal), the desired value for the total mass of the system (solid line) and the time evolution of (M(m)), the current total mass (in dashed line). The total mass convergence presents a similar behaviour to the reported simulated cases. Fig. 5. Real NMB response (upper plot) and drug dose profile (lower plot) administered to the real patient ( ˆ Θidentified with the Bayesian algorithm from t0until t1) that result from the control with the Total Target Control Algorithm. The ’dot’ in the xx-axis indicates time instant t0computed by OLARD. Note that NMB signal is corrupted with sensor noise. IV. CONCLUSIONS This paper presents an integrated algorithm for automatic drug dose administration referred as Total Mass Target Controlled Infusion (TCI) driven by parametric estimation and a robust mass control of compartmental systems with parameters uncertainties. The performance of the developed strategy was evaluated for the NMB (neuromuscular ) control. This particular appli- Fig. 6. Parameter evolution from minute 10 until t1= 75 min for the simulation control essay which results are represented in Fig. 5. Refinements of the parameterization were carried out from t0= 40 min until t1. The values are normalized with the population parameter values of models present in M. Fig. 7. Total mass evolution for the real control case. In solid line is represented the target level in total mass (M∗(ˆ Θfinal)); in dashed line is represented the time evolution of the system total mass (M(m)); in dots is represented the different values of M∗(ˆ Θk)in the transient phase. cation enables the true evaluation of the control system since it has a reliable effect sensor and is widely accepted in clinical environment. The proposed Total Mass TCI strategy comprises three different stages. After the administration of the initial bolus the dose to be administrated is zero until the recovery of the signal is detected. First, this identification was carried out by a developed empirical algorithm, OLARD (On Line tuned Algorithm for Recovery Detection). Secondly, during the transient phase, i.e., until the steady state of the reference attains the target value, an identification method together with an Input Variance Control (IVC) iterative implemented algorithm is used to determine the optimal dose profile. The identification method is a Bayesian algorithm based in sparse measurements of the signal and proved to obtain reliable estimation for the model parameters. Finally, driven by the parameter estimation a Total Mass Control (TMC) scheme for compartmental systems which proved to have robustness properties under the presence of uncertainties is applied to maintain the level in a constant target value. The whole strategy was applied to a bank of models, and demonstrated to be adequate to the design of individualized dose regimen both on open and closed loop systems, even in the presence of noise and uncertainties. The main contributions of this paper, namely the developed identification and control algorithms and the results obtained encourage a stake on dedicated improvements in order to obtain an enhanced and robust version of a Target Control Infusion prototype. ACKNOWLEDGMENT This work was supported by Fundac¸˜ ao para a Ciˆ encia e Tecnologia (FCT) through Project IDEA (reference PTDC/EEAACR/69288/2006). The work of R. Sebasti˜ ao is supported by FCT under the PhD Grant SFRH/BD/41569/2007. REFERENCES [1] P. J. Lago, “Open-loop stochastic control of pharmacokinetic systems: a new method for design of dosing regimens,” Comp. and Biom. Research, vol. 25, pp. 85–100, 1992. [2] D. S. Bayard, M. H. Milman, and A. Schumitzky, “Design of dosage regimens: a multiple model stochastic control approach,” Int. Journal of Bio-Medical Computing, vol. 36, pp. 103–115, 1994. [3] A. R. Absalom and M. Struys, Overview of Target Controlled Infusions and Total Intravenous Anaesthesia. Gent: Academia Press, 2007. [4] C. H. Ting, R. H. Arnott, D. A. Linkens, and A. Angel, “Migrating from target-controlled infusion to closed-loop control in general anaesthesia,” Computer Methods and Programs in Biomedicine, vol. 75, pp. 127–139, 2004. [5] R. Sebasti˜ ao, M. M. Silva, S. Esteves, J. Gama, and T. Mendonc¸a, “Online Tuned Algorithm for Recovery Detection (OLARD),” accepted for presentation in the 2nd World Congress of Total Intravenous Anaesthesia - TCI, (Berlin, Germany), 2009. [6] T. Mendonc¸a and P. Lago, “PID control strategies for the automatic control of neuromuscular blockade,” Control Engineering Practice, vol. 6, no. 10, pp. 1225–1231, 1998. [7] P. Lago, T. Mendonc¸a, and L. Gonc¸alves, “On-line autocalibration of a PID controller of neuromuscular blockade,” in IEEE Int. Conference on Control Applications, (Trieste, Italy), pp. 363–367, 1998. [8] T. C. Jannett and S. Aragula, “Simulation of Adaptive Control of Theophylline Concentration,” IEEE Control Systems, pp. 32–37, 1992. [9] H. Alonso, T. Mendonc¸a, and P. Rocha, “A hybrid method for parameter estimation and its application to biomedical systems,” Computer Methods and Programs in Biomedicine, vol. 89, pp. 112–122, 2008. [10] W. J. Rugh, Linear System Theory. USA: Prentice Hall Information and System Sciences Series, 1996. [11] H. T. Toivonen, “Variance constrained self-tuning control,” Automatica, vol. 19(4), pp. 415–418, 1983. [12] C. Sousa, T. Mendonc¸a and P. Rocha,“Total mass control in uncertain compartmental systems,” Proceedings of the 8th Portuguese Conference on Automatic Control, (Vila Real, Portugal), 2008. [13] T. Mendonc¸a, H. Magalh˜ aes, P. Lago, S. Esteves, “Hippocrates: a robust system for the control of neuromuscular blockade,” Journal of Clinical Monitoring and Computing,vol. 18, pp. 265–273, 2004.