c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

journal homepage: www.intl.elsevierhealth.com/journals/cmpb

Local identifiability and sensitivity analysis of neuromuscular blockade and depth of hypnosis models M.M. Silva a,b,c,∗ , J.M. Lemos d , A. Coito d , B.A. Costa d , T. Wigren b , T. Mendonc¸a a,c a

Departamento de Matemática, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal b Division of Systems and Control, Department of Information Technology, Uppsala University, Box 337, SE-751 05 Uppsala, Sweden c Center for Research and Development in Mathematics and Applications (CIDMA), Universidade de Aveiro, Campus Universitário de Santiago, 3810-193 Aveiro, Portugal d INESC-ID, Instituto Superior Técnico, Technical University of Lisbon, Rua Alves Redol 9, 1000-029 Lisboa, Portugal

a r t i c l e

i n f o

a b s t r a c t

Article history:

This paper addresses the local identifiability and sensitivity properties of two classes of

Received 5 July 2012

Wiener models for the neuromuscular blockade and depth of hypnosis, when drug dose pro-

Received in revised form

files like the ones commonly administered in the clinical practice are used as model inputs.

22 June 2013

The local parameter identifiability was assessed based on the singular value decomposi-

Accepted 20 July 2013

tion of the normalized sensitivity matrix. For the given input signal excitation, the results show an over-parameterization of the standard pharmacokinetic/pharmacodynamic mod-

Keywords:

els. The same identifiability assessment was performed on recently proposed minimally

Anesthesia

parameterized parsimonious models for both the neuromuscular blockade and the depth

Identifiability

of hypnosis. The results show that the majority of the model parameters are identifiable

Minimally parameterized models

from the available input–output data. This indicates that any identification strategy based

PK/PD models

on the minimally parameterized parsimonious Wiener models for the neuromuscular block-

Sensitivity

ade and for the depth of hypnosis is likely to be more successful than if standard models

Wiener models

are used. © 2013 Elsevier Ireland Ltd. All rights reserved.

Abbreviations: BIS, bispectral index; DoH, depth of hypnosis; EEG, electroencephalogram; MPP, minimally parameterized parsimonious; NMB, neuromuscular blockade; PD, pharmacodynamics; PEM, prediction error method; PK, pharmacokinetics; TOF, train-of-four. ∗ Corresponding author at: Division of Systems and Control, Department of Information Technology, Uppsala University, Box 337, SE-751 05 Uppsala, Sweden. Tel.: +46 184712849. E-mail addresses: [email protected], [email protected] (M.M. Silva), [email protected] (J.M. Lemos), [email protected] (A. Coito), [email protected] (B.A. Costa), [email protected] (T. Wigren), [email protected] (T. Mendonc¸a). 0169-2607/$ – see front matter © 2013 Elsevier Ireland Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cmpb.2013.07.020

24

1.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

Introduction

This paper addresses the local identifiability and sensitivity properties of two classes of Wiener models for the neuromuscular blockade and the depth of hypnosis of patients subject to general anesthesia.

1.1.

General anesthesia

General anesthesia is a reversible drug-induced state in which the patient’s muscle relaxation, analgesia and hypnosis are guaranteed [1]. The degree of muscle relaxation can be measured from the electromyography evoked at the hand of the patient by electrical stimulation of the adductor pollicis muscle to supramaximal train-of-four (TOF) stimulation of the ulnar nerve. In this paper, the measurement of the neuromuscular blockade (NMB) was chosen as the first single response from a TOF stimulation calibrated by a reference twitch [2]. Consequently, the quantification of the NMB is normalized between 0% (total paralysis) and 100% (full muscular activity). Analgesics aim to provide pain relief. At present, a quantitative and reliable index for the direct measurement of pain in patients has not yet been widely accepted and validated. Hypnotics control the level of unconsciousness and must guarantee absence of recall of intra operative events after surgery. Since most of the hypnotics and analgesics interact, the effect of the administration analgesics in the depth of hypnosis (DoH) modeling is usually taken into account. Several univariate parameters computed using the raw data from the electroencephalogram (EEG) are used to assess the DoH. The bispectral index (BIS) [3] is the most commonly used index by anesthesiologists and researchers in the field to infer the patients’ DoH. The BIS ranges from 0, equivalent to a flat line EEG, to 100 (or 97.7 depending on the software version implemented in the BIS monitor), equivalent to a fully awake state. After the induction of a general anesthesia, the BIS should lay between 40 and 65 [4].

1.2.

Background

From an engineering point of view, a patient under anesthesia may be envisaged as a black box from which the only accessible measurements are the amount of drugs that are administered (the system inputs), and the physiological responses that are monitored (the system outputs). It is too seldom that intermediate variables such as drug concentrations in different parts of the body are measured. Consequently, in order to predict the physiological responses or unobserved concentrations that result from different inputs or external perturbations [5], mathematical models of the effect of drugs have to be built. A major bottleneck of the accuracy of these predictions is the estimation of the model parameters from the measured signals. The possibility to uniquely recover the unknown parameters from input–output data is denoted identifiability. A priori global identifiability of linear systems has been extensively studied [6]. However, the systems describing the effect of drugs in the human body are usually nonlinear. Therefore, issues like identifiability of those

models and the convergence of the estimates have few general answers [7]. Even though being important for display, warning, decision analysis and outcome comparison [8], a significant use of the models describing the effect of drugs in the human body is to design closed-loop controllers. When identifying systems operating in closed-loop, two important aspects need to be jointly considered [9]: the excitatory properties of the controlled input signals and the model structure. This follows since it is known that a more complicated model with more parameters puts higher requirements on the inputs than do simple models with few parameters [10]. In the case of drug delivery in anesthesia, the drug dose profiles (input signals) have to follow standard administration protocols, and cannot be arbitrarily selected for best identification experiments or probing excitation [8]. In intravenous anesthesia the common practice is to administer a bolus dose (that approximates an impulse from the mathematical point of view) at the beginning of the surgery, followed by a sequence of steps of different amplitudes depending on the patients’ needs. Given this poorly exciting profile of the (open-loop and closed-loop) input signals in this application, the choice of the model structure is crucial to guarantee model identifiability. By looking at the properties of the normalized sensitivity matrix [11], the sensitivity and identifiability properties of the model at a given point in the parameter space can be studied together with the excitatory properties of the input signals. Hence, what is here at stake is the combined effect of the input signal excitation and the number of parameters to be identified in the models. This identifiability issue was implicitly mentioned in [12], where the authors acknowledge that the main identification challenge lies in that the input–output data is not informative enough to robustly distinguish between two parameters in the standard pharmacokinetic/pharmacodynamic (PK/PD) model used to describe the effect of the hypnotic propofol in anesthesia. Some studies have been recently carried out to explicitly assess the identifiability properties of the standard PK/PD models describing the effect of drugs in anesthesia. Regarding the NMB, the results of a model sensitivity analysis in [13] show that half of the parameters in the standard PK/PD model for the muscle relaxant atracurium are potentially unidentifiable. In the paper, this information is useful to modify the weights of a Bayesian cost function for the identification of the model parameters. This approach was further extended in [14] to the standard DoH model, using data from sedation cases. The conclusions of the local parameter identifiability assessment in [14] are that the standard PK model is not fully identifiable, which suggests the use of other, simpler, models. Minimally parameterized parsimonious (MPP) models for these two components of anesthesia were recently proposed [15,16], and no identifiability study on those models has been carried out yet.

1.3.

Contributions and paper structure

The first contribution of this paper is to investigate the sensitivity with respect to the model parameters of both the standard PK/PD and the MPP Wiener models for the NMB and

25

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

the DoH. This analysis provides information about the time periods in which data is more informative for parameter estimation, and is the basis for the local identifiability analysis that is also performed in this paper for both the standard PK/PD and the MPP models, this being the second contribution of the paper. A central concept in this framework is the fact that the identifiability properties of the models is a structural problem about the models under consideration and, most importantly, is independent on which method is chosen to estimate the model parameters. The remainder of this paper is organized as follows. Section 2.1 presents both the standard PK/PD and the MPP models for the NMB. Section 2.2 presents both the standard PK/PD and the MPP models for the DoH. The methods used to assess sensitivity and identifiability of the four models are described in Section 3. Section 4 presents the results, and Section 5 draws the conclusions.

A possible state-space realization of (1) is

⎡

x˙ 1 (t)

⎤

⎡

−1

⎢ ⎥ ⎢ ⎢ x˙ 2 (t) ⎥ ⎢ 0 ⎢ ⎥=⎢ ⎢ ⎥ ⎢ ⎣ c˙ (t) ⎦ ⎣  c˙ e (t)

0

0

0

0

−2

0

0



−

0

1/

⎤⎡

x1 (t)

⎤ ⎡

a1

⎤

⎥⎢ ⎥ ⎢ ⎥ ⎥ ⎢ x2 (t) ⎥ ⎢ a2 ⎥ ⎥⎢ ⎥ + ⎢ ⎥ ua (t), ⎥⎢ ⎥ ⎢ ⎥ ⎣ ⎦ c(t) ⎦ ⎣ 0 ⎦ 0 −1/

0

ce (t)

(2) where xi (t), i = 1, 2 are nonnegative state variables such that the atracurium plasma concentration cp (t) equals x1 (t) + x2 (t). The state variable c(t) is used as an intermediate variable. Note that in the clinical practice none of the intermediate signals c(t), cp (t) and ce (t) are measured. The nonlinearity in the Wiener structure relates the effect concentration ce (t) to the effect of the drug as quantified by the measured NMB y(t) [%] and is usually modeled by the Hill function [18] as 

2.

y(t) =

Models in anesthesia

A common feature among the majority of the models that are available in the literature to describe the effect of drugs in the human body is that they are usually constituted by a linear block followed by a nonlinear function. The linear block models the distribution of the drug in the body i.e. the relation between the administered drug and its plasma and effect concentration, while the nonlinearity models the relationship between the drug effect concentration and its measured effect. Regardless of the choice of the specific model structure inside each block, a linear dynamic model followed by a static nonlinearity is called a Wiener model [17]. The following subsections present the specific models for both the linear dynamics and static nonlinearity in the Wiener models for the NMB and the DoH.

100C50



Neuromuscular blockade

2.1.1.

NMB: standard PK/PD model

The standard PK/PD model for the effect of the nondepolarizing muscle relaxant atracurium in the NMB is presented in this subsection. The NMB input–output relationship is commonly described as a Wiener model [18,19], and consists of a linear dynamic part that can be seen as a compartmental model [20], followed by a static nonlinearity. In the frequency domain, the linear dynamics using macro constants may be modeled in continuous-time as

Ce (s) =

 a 1 s + 1

+

a2 s + 2

 

1/ Ua (s), s +  s + 1/

(1)

where Ce (s) is the Laplace transform of the effect concentration of atracurium ce (t) [␮g ml−1 ], Ua (s) is the Laplace transform of the atracurium infusion rate ua (t) [␮g kg−1 min−1 ] (system input), and {ai [kg ml−1 ], i [min−1 ]}i=1,2 , [min−1 ] and [min] are patient-dependent parameters.

(3)

,

where C50 [␮g ml−1 ] and  (dimensionless) are also patientdependent parameters. The total number of parameters in the NMB standard PK/PD model characterizing each individual patient response is hence eight (see (27) for a full list).

2.1.2. model

NMB: minimally parameterized parsimonious

Mainly due to the high number of parameters to be identified in the standard PK/PD model, a new MPP model for the NMB was proposed in [15]. In the frequency domain, the linear dynamics of the NMB MPP model may be described in continuous-time by Ya (s) =

2.1.



C50 + ce (t)

˛ k1 ˛ k2 ˛ Ua (s), s + ˛ s + k1 ˛ s + k2 ˛

(4)

where Ya (s) is the Laplace transform of the continuous-time output ya (t) of the model linear dynamic block and Ua (s) is the Laplace transform of the atracurium infusion rate ua (t) [␮g kg−1 min−1 ]. Here the patient-dependent parameter is ˛. For the muscle relaxants atracurium and rocuronium, the parameters ki , i = 1, 2 are chosen to be constant and equal to 4, and 10, respectively, as explained in [21,22]. Further studies are needed to determine these constants to other muscle relaxants. A state-space (compartmental) realization of (4) is

⎡

x˙ 1 (t)

⎤

⎡

−k2 ˛

⎢ ⎥ ⎢ ⎣ x˙ 2 (t) ⎦ = ⎣ k1 ˛ y˙ a (t)

0

0

0

⎤⎡

−k1 ˛

0

⎥⎢ ⎥ ⎢ ⎥ ⎦ ⎣ x2 (t) ⎦ + ⎣ 0 ⎦ ua (t).

˛

−˛

x1 (t)

ya (t)

⎤ ⎡

k2 ˛

⎤ (5)

0

The major advantage of this new model for the linear dynamic block is that it reduces the inter-patient variability to only one parameter, in contrast with six parameters in the standard compartmental model for the NMB. As described in Section 2.2, the same happens with the modeling of the DoH.

26

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

The static nonlinearity is the same as in the standard model 

y(t) =

100C50

 C50



+ ya (t)

(6)

,

where  is also a patient-dependent parameter that adapts the shape or static incremental gain of (6). The parameter C50 , different from the one in (3), is chosen to be constant as explained in [21]. The total number of parameters in the NMB MPP model characterizing each individual patient response is hence two (see (32)).

2.2.

DoH: standard PK/PD model

The standard PK/PD model to describe the joint effect of the hypnotic propofol and the analgesic remifentanil in the DoH is described in this subsection. Both for propofol and remifentanil, a three-compartment mammilary model is normally used to explain the linear distribution of each drug in the different theoretical compartments of the human body [24,25]. By direct deduction from mass balances between compartments, and assuming that for each compartment i at time t, an amount xi (t) [mg] of drug is present, the state space representation that is commonly used becomes

⎡

x˙ 1 (t)

=

Up (t) , Up (t) + Ur (t)

(8)

where, by definition,  ranges from 0 (remifentanil only) to 1 (propofol only). In the remaining of this paper the subscript “p” refers to propofol and the subscript “r” to remifentanil. To calculate (8) both effect concentrations cep (t) and cer (t) are first normalized with respect to their concentration at half the maximal effect (C50p and C50r , respectively) as

Depth of hypnosis

Anesthetic drugs are often combined because they interact synergistically to create the anesthetized state [23]. Due to this fact, during a clinical procedure, the desired DoH is often maintained by the combined administration of a hypnotic and an analgesic. The main common feature between the NMB and the DoH is that both are usually described by Wiener models.

2.2.1.

The nonlinear supra-additive model proposed in [23] is the most commonly used model to describe the joint effect of propofol and remifentanil in the DoH. The potency of the drug mixture is hence modeled as

⎤

⎡

⎢ ⎥ ⎢ ⎢ x˙ 2 (t) ⎥ ⎢ ⎢ ⎥ ⎢ ⎢˙ ⎥=⎢ ⎣ x3 (t) ⎦ ⎣

−(k10 + k12 + k13 )

k21

k31

0

k12

−k21

0

0

k13

0

−k31

c˙ e (t)

ke0 /V1

0

0

⎤⎡

x1 (t)

⎤

⎥⎢ ⎥ ⎥ ⎢ x2 (t) ⎥ ⎥⎢ ⎥ ⎥⎢ ⎥ 0 ⎦ ⎣ x3 (t) ⎦ −ke0

ce (t)

⎡ ⎤ 1

⎢0⎥ + ⎢ ⎥ u(t), ⎣0⎦

(7)

0 where u(t) [mg ml−1 min−1 ] is the drug infusion rate (either propofol or remifentanil), V1 [ml] is the volume of the central compartment (compartment 1), k10 [min−1 ] is the clearance of the drug modeled by an outer flow of the compartment 1, and kij [min−1 ] are transfer coefficients from compartment i to compartment j. The dynamics of the effect concentration of the drug ce (t) [mg ml−1 ] is given by the last row in (7) [26], where ke0 describes the rate of removal of drug from the effectsite out of the body. In clinical practice, neither xi (t) nor ce (t) can be measured.

Up (t) =

cep (t) , C50p

Ur (t) =

cer (t) . C50r

(9)

The nonlinear concentration–response relationship for any ratio of these two drugs can then be described by a generalized Hill function as y(t) =

y0 , 1 + ((Up (t) + Ur (t))/U50 ())

(10)

where y0 is the effect at zero concentration,  controls the steepness of the nonlinear concentration–response relation, and U50 () is the normalized concentration associated with 50% of the maximum effect of both drugs at ratio . In [23] the quadratic polynomial U50 () = 1 − ˇ + ˇ2

(11)

was proposed for the expression of U50 (). In common clinical practice, the parameters , ˇ, y0 , C50p and C50r in (9)–(11), and the linear parameters kij , ke0 and V1 in (7) are usually estimated based on population model distributions and on the clinicians’ expertise. The total number of patient-dependent parameters in the DoH standard model is hence nineteen (see (29) for a full list).

2.2.2.

DoH: minimally parameterized parsimonious model

The MPP Wiener model describing the joint effect of the hypnotic propofol and the analgesic remifentanil in the DoH is described in this subsection. In [16], a third-order continuoustime model was proposed for the linear dynamics of both propofol and remifentanil. In the frequency domain, propofol linear dynamics may be modeled in continuous-time by Yp (s) =

d1  d2   Up (s), s +  s + d1  s + d2 

(12)

where Yp (s) is the Laplace transform of the continuous-time output yp (t) from the linear dynamic block of the model for propofol; Up (s) is the Laplace transform of the propofol infusion rate up (t) (input signal);  is a patient-dependent parameter; and the parameters di , {i = 1, 2} are chosen to be constant and equal to 9 and 10, respectively, as explained in [27].

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

In the frequency domain, remifentanil linear dynamics may be modeled in continuous-time by Yr (s) =

 l1  l2  Ur (s), s +  s + l1  s + l2 

(13)

27

expansion of y(t, ) with respect to around 0 and neglecting high order terms, the increments ı and ıy(t) are related by: ıy(t) = Sy (t)ı ,

(17)

where where Yr (s) is the Laplace transform of the continuous-time output yr (t) of the linear dynamic block of the model for remifentanil; Ur (s) is the Laplace transform of the remifentanil infusion rate ur (t) (input signal);  is a patient-dependent parameter; and the parameters li , {i = 1, 2} are chosen to be constant and equal to 2 and 3, respectively, as explained in [27]. Both linear dynamics (12) and (13) may be represented in a compartmental state-space realization as (5). Using as starting point the standard model (10), a new formulation for the nonlinearity was proposed in [16], being y(t) =

y0 1 + ((yr (t)/C50r ) + m(yp (t)/C50p ))



,

(14)

where m and  are patient-dependent parameters and y0 is equal to 97.7. C50p and C50r are propofol and remifentanil normalizing constants (different from the ones in (9)), respectively, [16]. From this, the total number of patient-dependent parameters in the DoH MPP model is eleven (see (33)).

3.

Methods

3.1.

Sensitivity analysis and local identifiability



∂y ∂y ··· ∂ 1 ∂ np



(15)

Sy,j (t) =

∂y(t, ) , ∂ j

(19)

and differentiating (16) with respect to a parameter j , the equation of the sensitivity of the output to a parameter j becomes ∂ ∂h ∂x ∂h y(t, ) = + , ∂ j ∂x ∂ j ∂ j

(20)

where Sx,j (t) = ∂x/∂ j , a column vector with the dimension n of the state space vector, is the state sensitivity with respect to a generic parameter j . By differentiating (15), interchanging the order of derivatives and using the chain rule, it follows that (21)

From the sensitivity computation, a local identifiability analysis is performed around the nominal value of the parameters. For that purpose, let the output samples of y be stacked into a single column Y and let ıY be a perturbation in the model output induced by a small parameter perturbation ı . These quantities are related by ıY = Sı ,

(22)

where each row of the sensitivity matrix S corresponds to Sy (t), as computed in (18), for a given time instant t. A further step is taken in order to eliminate the parameter scaling effects. This is done by defining the normalized parameter perturbation 

−1

ı =[diag( )]

ı .

(23)

With this definition it follows that ıY = Sı ,

and a nonlinear output equation written as

(18)

is the sensitivity row vector of dimension equal to the number of parameters np . This follows since y is a scalar in this paper. As mentioned before, the components of the sensitivity vector depend on time. From the definition of sensitivity,

∂A( ) ∂B( ) d x(t) + A( )Sy,j (t) + u(t). S (t) = ∂ j ∂ j dt x,j

Sensitivity analysis assesses how the variation of a system parameter affects the system output. The sensitivity of the output is time-dependent and is given by a function that satisfies a differential equation. This allows to find when in time a certain parameter has its greatest effect. If a small change in a parameter results in relatively large changes in the output, the output is said to be sensitive to that parameter. If parameters have a small sensitivity, they tend to be hard to estimate since a small deviation does not induce an evident change in the output response of the system, and may be masked by noise. This subsection presents the methods used for the sensitivity analysis and local identifiability assessment of the NMB and DoH models. As discussed in Sections 2.1 and 2.2 the NMB and the DoH input–output relationships can be modeled by a continuoustime linear dynamic part, which is represented by ˙ = A( )x(t) + B( )u(t), x(t)

∂y Sy (t) = = ∂

(24)

where the normalized sensitivity matrix S is given by y(t, ) = h(x(t), ).

(16) S = Sdiag( ).

Let 0 be the vector of nominal parameters and the perturbed one. The perturbation ı = − 0 induces a perturbation ıy(t) = y(t, ) − y(t, 0 ) on the output. Performing a Taylor

(25)

The advantage of combining all the output samples into a single matrix is the possibility to verify if there are linear

28

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

dependencies between parameters that make the parameters not (locally) identifiable. The rank of the matrix provides information about the effect of parameter changes in the output. If the sensitivity matrix does not have full rank, then some parameters are unidentifiable because different nonzero combinations of their increments yield the same output increment. In practice, this is evaluated by performing a singular value decomposition (SVD) of the sensitivity matrix S (for absolute parameter variations) or of the normalized sensitivity matrix S (for relative parameter variations) and looking at its singular values. The SVD of an m × n matrix M is a factorization of the form M = U˙V,

The inter-individual variabilities of the parameters C50p , C50r ,  and ˇ were modeled using a log-normal variance model with mean and standard errors equal to 3.40 ± 0.75 [g ml−1 ], 8.91 ± 2.35 [ng ml−1 ], 4.29 ± 0.98, and 1.69 ± 0.42, respectively, as in [32]. The mean values and standard deviations for the nominal parameters of the DoH standard PK/PD model are shown in Table A.2.

3.2.2.

Minimally parameterized parsimonious models

Since no general database of parameters for the NMB and DoH MPP models has been published until the writing of this paper, identification methods have to be used to obtain nominal values of the parameters in those models. The prediction error method (PEM) proposed in [21] was used. In general, the PEM determines the parameter vector so that the prediction error

(26)

where U is a m × m matrix, ˙ a m × n diagonal matrix with nonnegative real numbers (the singular values of M) on the diagonal and V a n × n matrix. The identifiability of the parameters present in the model may then be obtained from the singular vectors of the sensitivity matrix. The entries of the singular vector corresponding to the lower singular values are high for parameters that may be poorly identifiable. These parameters can e.g. be fixed a priori, reducing the complexity of the identification task and eliminating useless degrees of freedom [28].

ε(t, ) = y(t) − yˆ (t, )

(29)

becomes as small as possible. In (29), y(t) is the measured output (NMB or DoH) and yˆ (t, ) is the predicted output based on the model that depends on the parameter vector . The criterion to be minimized in the search, [10], is 1 2 ε (t, ), N N

3.2.

V( ) =

Nominal parameters

(30)

t=1

As mentioned in Section 3.1, in order to compute the sensitivity and to perform a local identifiability analysis, the nominal values of the model parameters are needed.

3.2.1.

where N is the total number of data points. The minimization of (30) was performed using the numerical Gauss–Newton algorithm, [10],

Standard PK/PD models

ˆ (k+1) = ˆ (k) − ˇ(V  ( ˆ (k) ))

For the NMB, taking into account the parameterization in Section 2.1.1, and in order to cover a wide range of behaviors, a bank M of models Mi = { i , i = 1, . . ., 100} was generated by using the probabilistic model discussed in [29]. The parameter vector for the NMB standard model is given by = [a1

1

a2

2

 C50



]T .

k12p

k13p

k21p

k31p

ke0p

(27)

V1p

k10r

k12r

k13r

A bank N of models Ni = { i , i = 1, . . ., 50} was generated for patients 18–70 years old, 160–180 cm height, 60-80 kg weight, 25 males, 25 females, according to the regressions in [31, pp. 89–90]. As it also happens with the MPP models, some of the parameters in the PK/PD model for the DoH are fixed (see Table A.2).

T

ˆ (k) ) , V  (

(31)

ˆ (k) denotes the kth iteration point in the search and where ˇ is a diagonal matrix used to control the step length. The derivatives of V( ) can be found in [21]. As mentioned in Section 2.1.2, ˛ and  are the patientdependent parameters in the NMB MPP model. Since the remaining constants are common to all patients, they are not considered in the identification step. Hence, for the NMB, the PEM was run in a database of 60 real clinical cases collected in the surgery room [21] for = [˛ ] in (29)–(31), and the nominal values of the parameters ˛ and  were obtained, forming the new bank P of models Pi = { i , i = 1, . . ., 60} of nominal values for the parameters in the NMB MPP model. The parameter vector for the sensitivity and identifiability assessment of the NMB MPP model stands as

The mean values and standard deviations for the nominal parameters of the NMB standard PK/PD model are shown in Table A.1. For the DoH, in the clinical practice, population models are applied to infer the DoH response of each patient to a certain drug administration profile. For propofol, the parameters from Schnider [30] are usually used, whereas for remifentanil, the parameters from Minto [25] are commonly applied. The joint parameter vector for the DoH standard model accounting for the effect of both propofol and remifentanil stands as = [k10p

−1

= [k1

k21r

k31r

k2

ke0r

˛ 

V1r

C50 ]T .

C50p

C50r

(32)



ˇ

y0 ]T .

(28)

The mean values and standard deviations that were obtained for the nominal parameters of the NMB minimally parameterized parsimonious model are shown in Table A.3. Similarly, as mentioned in Section 2.2.2, , , m and  are the four patient-dependent parameters in the DoH MPP model. Since the remaining constants are common to all patients,

29

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

1500 ua (t)[ µg/k g/min]

they are not considered in the identification step. Hence, for the DoH, the PEM was run in a database of 23 real cases collected in the surgery room [21] and the nominal values of the parameters vector = [  m ]T were obtained via (29)–(31), and used in the new bank Q of models Qi = { i , i = 1, . . ., 23} of nominal values for the parameters in the DoH MPP model. The parameter vector for the sensitivity and identifiability assessment of the DoH MPP model stands as

... 20 15 10 5 0 0

50

100 150 Time [min]

200

240

Fig. 1 – Input (atracurium rate) for the NMB simulations. d2

l1

l2



 C50p

C50r

m 

y0 ] .

(33)

The mean values and standard deviations that were obtained for the nominal parameters of the DoH minimally parameterized parsimonious model are shown in Table A.4. For consistency purposes it is important to assess sensitivity and identifiability properties of both the patient-dependent parameters and the ones that were chosen to be constant among patients. The results of such assessment may either support the choices that were made in [15,16], or point to improvements or changes regarding the parameters that have been considered to span the patient dynamic variability of behaviors and that were consequently identified in the referred works. Due to the fact that both the PK/PD models and the MPP models in [18,30,25,23,15,16], respectively, were derived using data collected from adults under general anesthesia, the sensitivity and local identifiability analysis presented in this paper is only valid for a population of adults. As a result, the conclusions taken from the work described in this manuscript cannot be directly extrapolated to elderly or children until further results are published.

4.

Results

The methodology described in Section 3.1 was applied using each set of nominal parameters in the four banks of models described in Section 3.2. Hence, for each bank of models, a normalized sensitivity matrix S(i) was computed for each case i. An

100 75 y(t)[%]

= [d1

T

50 25 0 0

50

100 150 Time[min]

200

240

Fig. 2 – Average of the outputs of the 100 standard PK/PD models for the NMB in the bank M, using the input of Fig. 1.

sensitive with respect to the parameters correspond to the drop and the recovery after the initial bolus (see Fig. 2). This was expected since the transient periods reveal the system dynamics more than stationary ones. Fig. 3(b) shows that the parameters corresponding to the largest output sensitivity are a2 , , 2 , and C50 . During the recovery from the initial bolus (between minute 50 and 100), the model is more sensitive to  and . On the opposite, the output is not so sensitive to parameters a1 , 1 and . Fig. 4 shows the singular values of the normalized sensitivity matrix S for the NMB standard PK/PD model. Since all the singular values are within a reasonable range, all parameters in the model are identifiable, although not all with equal

N

(a) 1000

a1 λ1 a2 λ2 λ C50 γ τ

0 Sy (t)

averaged sensitivity matrix S = (1/N) · i=1 S(i) , with N denoting the number of models in the bank, was also computed using the normalized sensitivity matrices S(i) . In order to perform the identifiability analysis, S was decomposed in singular values. The results obtained are presented in this section.

−1000 −2000 −3000

NMB results

This subsection presents the results for the NMB standard PK/PD and MPP models described in Section 2.1. For the sensitivity/identifiability analysis assessments, the same atracurium profile, shown in Fig. 1, was applied to all models in the banks.

−4000 0

50

100 150 Time [min]

200

240

a1 λ1 a2 λ2 λ C50 γ τ

(b) 100 50 Sy (t)

4.1.

0 −50

4.1.1.

Standard PK/PD model

Fig. 2 shows the time evolution of the average of the outputs of the 100 standard PK/PD models for the NMB in the bank M. Fig. 3(a) shows the time-varying sensitivity of the NMB standard PK/PD model output with respect to the parameter vector in (27). The periods in which the output is more

−100 0

25

50 Time [min]

75

100

Fig. 3 – Sensitivity results for the NMB standard PK/PD model. (a) Output sensitivity with respect to the model parameters. (b) Zoom of (a).

30

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

0

0

log10 (SV )

log10 (SV )

−1 −2 −3

−1

−2

−4 3

5 SV Index

7

−3 1

8

Fig. 4 – Singular values of S for the NMB standard PK/PD model.

100

y(t)[%]

75 50 25 0 0

50

100 150 Time [min]

200

Sy (t)

k1 k2 α γ C50

−1000 −2000 50

100 150 Time [min]

200

240

(b) 20

k1 k2 α γ C50

Sy (t)

10 0 −10 −20 0

50

100 150 Time [min]

200

60 ...

4

5

10 5 0 0

ur (t)[mg/min ]

significance. In order to reduce the complexity of a possible identification task, the lowest singular values, namely the seventh and eighth eigenvalues, can be considered negligible, eliminating useless degrees of freedom. By looking at the singular vectors of the sensitivity matrix, shown in Table B.5, the entries of the seventh and eighth singular vectors are high (values in bold) for parameters that may not be easily identifiable: in this case a1 and a2 .

0

3 SV Index

20

40 60 Time [min]

80

100

20

40 60 Time [min]

80

100

(b) 0.12

Fig. 5 – Average of the outputs of the 60 minimally parameterized parsimonious models for the NMB in the bank N, using the input of Fig. 1.

−3000 0

(a)

240

(a) 1000

2

Fig. 7 – Singular values of S for the NMB minimally parameterized parsimonious model.

up (t)[mg/min ]

−5 1

240

Fig. 6 – Sensitivity results for the NMB minimally parameterized parsimonious model. (a) Output sensitivity with respect to the model parameters. (b) Zoom of (a).

... 0.01 0.005 0 0

Fig. 8 – Inputs for the DoH simulations. (a) Input (propofol infusion rate). (b) Input (remifentanil infusion rate).

4.1.2.

Minimally parameterized parsimonious model

As shown in Fig. 6(a), the output (as shown in Fig. 5) of the NMB MPP model is very sensitive to the parameter ˛. The periods during which major changes in the parameters sensitivity occur correspond to the drop and recovery after the initial bolus (Fig. 6(b)). The parameter that corresponds to the largest sensitivity is ˛, even though C50 and  are also quite sensitive. On the contrary, the output is little sensitive to the parameters k1 and k2 . Fig. 7 shows the singular values of the normalized sensitivity matrix S for the NMB MPP model. From the dispersion of the singular values, the fifth eigenvalue can be negligible. Similarly as before, by looking at the singular vectors of S, shown in Table B.6, the entries of the fifth singular vector are high (values in bold) for parameters that may not be identifiable: in this case k1 and k2 . These results support the choice of reducing the number of parameters in the modeling of the NMB input–output relationship as described in [15]. The fact that not all parameters are easily identifiable in the standard PK/PD model was considered as the driving argument in [15] to propose a MPP model able to cope with the limited amount of available data and

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

For the sensitivity/identifiability assessment the same propofol and remifentanil profiles, shown in Fig. 8, were applied to all models in the banks.

100

y(t)

80 60 40

4.2.1.

20

Fig. 10(a) and (b) shows the time-varying sensitivity of the DoH standard PK/PD model output (see Fig. 9) with respect to the parameter vector in (27). The parameters that correspond to the largest sensitivity are k31p , k21p , k10p , and k13p . On the opposite, after the induction phase, the output is not sensitive to the parameters k13r , V1p , V1r , C50p , C50r , , ˇ and y0 . The sensitivity of the output with respect to V1r , C50p , C50r , , ˇ is, however, high during the induction phase. Fig. 11 shows the singular values of S for the DoH standard PK/PD model. The fourteenth to nineteenth singular values are considered to have a negligible contribution to the rank of matrix S. By looking at the singular vectors of the sensitivity matrix, shown in Table B.7, the entries of the corresponding singular vectors are high (values in bold) for parameters that may not be identifiable: in this case k13p , k31p , k13r , k31r , ke0r , C50r , and ˇ. The same structural assessment was made on the Marsh model [33], using the corresponding [31, p. 89] set of parameters for the population mentioned in Section 3.2.1. If the fourteenth to nineteenth singular values are considered negligible, additionally to the parameters that may not be

0 0

20

40

60

80

100

Time [min]

Fig. 9 – Average of the outputs of the 50 standard PK/PD models for the DoH in the bank P, using the input of Fig. 8.

little excitatory pattern of the input. In this paper a quantitative evidence of the local identifiability properties of the NMB standard PK/PD model is presented. Moreover, in the MPP model, the less identifiable parameters are k1 and k2 which supports the choice of fixing them a priori in the MPP model as mentioned in Section 2.1.2.

4.2.

31

DoH results

This subsection presents the results for the DoH standard PK/PD and MPP models described in Section 2.2.

Standard PK/PD model

(a) 400 300 200 Sy (t)

100 0 −100 −200 −300 −400 0

20

40

Time [min]

60

80

100

(b) 60 40

Sy (t)

20 0 −20 −40 −60 0

5

10 Time [min]

15

20

Fig. 10 – Sensitivity results for the DoH standard PK/PD model. (a) Output sensitivity with respect to the model parameters. (b) Zoom of (a).

32

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

0

−2

log10 (SV )

log10 (SV )

0

−4 −6 1

3

5

7

9 11 13 SV Index

15

17

−2 −4

19

−6 1

Fig. 11 – Singular values of S for the DoH standard PK/PD model.

3

5 7 SV Index

9

11

Fig. 14 – Singular values of S for the DoH minimally parameterized parsimonious model.

100

y(t)

80 60 40 0

20

40 60 Time [min]

80

100

Fig. 12 – Average of the outputs of the 23 minimally parameterized parsimonious models for the DoH in the bank Q, using the input of Fig. 8.

identifiable in the Schnider model, the SVD assessment show that k10p , V1p and k10r may also not be identifiable. These results, not shown in this manuscript, corroborate and strengthen the findings for the standard PK/PD model for the DoH.

4.2.2.

Minimally parameterized parsimonious model

As shown in Fig. 13, the output (as shown in Fig. 12) of the DoH MPP model is most sensitive to the parameters , , C50r and . Fig. 14 shows the singular values of the matrix S for the DoH MPP model. From the dispersion of the singular values the two smallest can be considered negligible. By looking at the singular vectors of the sensitivity matrix, shown in Table B.8, the entries of the singular vector corresponding to the tenth and eleventh singular values are high (values in bold) for parameters that may not be identifiable: k1 , k2 , l1 , l2 , C50p , C50r and m.

Fig. 13 – Output sensitivity for the DoH minimally parameterized parsimonious model with respect to the model parameters.

This result is apparently not in accordance with the model reduction proposed in [16] where k1 , k2 , l1 , l2 , C50p and C50r were the parameters to be fixed in the DoH MPP model, specially in what concerns C50p and m. However, if a careful analysis on the Hill equation in (14) is done, this result is not surprising. Indeed, since it is only the quotient between m and C50p that affects yp (t), one parameter has to be fixed so that the other is identifiable. This compromise, and the fact that y0 may be considered as an identifiable parameter, should be taken into account for a further improvement of the modeling presented in [16].

5.

Conclusions

In this work a sensitivity and local identifiability assessment of the standard PK/PD models for the neuromuscular blockade (NMB) [20,18] and for the depth of hypnosis (DoH) [23–25] was presented. The results from the local identifiability analysis suggest an over-parameterization of the standard PK/PD models for both the NMB and DoH, when drug dose profiles similar to the ones commonly administered in the clinical practice are used as the model input. With regard to model-depth, the demand “as simple as possible and as complex as necessary” plays one important role here. If system identification algorithms are to be used in an anesthesia environment to tackle the intra-patients variability, a model where some parameters are difficult to uniquely be identified from the available input–output measurements is not adequate. The same methodology was used to assess the sensitivity and local identifiability properties of recently proposed MPP models for the NMB [15] and for the DoH [16]. The results obtained show that any identification strategy based on these models benefits from the use of few parameters since they are identifiable from the available input–output data. These results support the premise in [8] that, due to the fact that any specific model structure is especially superior, simple model structures that can capture the key characteristics are preferable.

33

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

Table A.1 – Nominal parameters of the NMB standard PK/PD model.

Mean Std

a1

1

a2

2



C50





0.01770 0.007247

0.3622 0.1254

0.005070 0.001206

0.03465 0.004547

0.9983 0.008438

0.6527 0.02519

4.179 0.5475

1.744 0.9917

Table A.2 – Nominal parameters of the DoH standard PK/PD model.

Mean Std

Mean Std

k10p

k12p

0.4383 0.0560

0.3381 0.0883

k13p

k21p

0.1960 0

0.680 0.0019

k31p

ke0p

V1p

k10r

k12r

0.0035 0

0.4560 0

4.27 0

0.5084 0.0248

0.3806 0.0776

k13r

k21r

k31r

ke0r

V1r

C50p

C50r



ˇ

y0

0.041 0.0029

0.2013 0.0258

0.0127 0.0033

0.5491 0.1103

4.8274 0.4059

3.1977 0.5938

0.0091 0.0032

4.2566 0.9854

0.0138 0.0028

97.7 0

Table A.3 – Nominal parameters of the NMB minimally parameterized parsimonious model. k1 Mean Std

k2

4 0

10 0

˛



0.04322 0.01048

2.790 0.9744

C50 3.244 0

Table A.4 – Nominal parameters of the DoH minimally parameterized parsimonious model. d1 Mean Std

9 0

d2 10 0

l1

l2

2 0

3 0

C50p 10 0

C50r 0.1 0

y0





m



97.7 0

0.08818 0.06763

0.5886 0.5245

2.862 1.855

1.400 1.075

Table B.5 – Singular values and corresponding singular vectors (10−3 ) for results in Fig. 4. SV index SV a1 1 a2 2  C50  

1 3925

2 733

3 86.0

4 29.3

5 17.3

6 6.42

7 0.221

8 0.0486

−0.0408 6.42 0.255 8.47 24.7 2.25 −999 37.8

−1.16 −17.6 8.38 279 955 92.9 24.9 −30.2

−1.12 −31.7 −6.97 −32.2 106 −993 0.383 12.7

−1.09 915 15.2 371 −90.3 −52.5 2.12 −121

−0.955 3.34 6.09 288 −54.5 −3.01 37.2 955

−12.8 402 −95.6 −832 255 45.3 11.9 264

−30.9 −24.7 − 995 90.2 −19.0 2.57 −0.935 −21.8

999 5.34 − 32.0 −6.92 3.76 −0.401 0.150 3.48

Mode of availability of software The software containing the sensitivity and local identifiability algorithms that were described in this paper is available on request by sending an email to Teresa Mendonc¸a ([email protected]).

Applications and Fundac¸ão para a Ciência e a Tecnologia, within project PEst-C/MAT/UI4106/2011 with COMPETE number FCOMP-01-0124-FEDER-022690, and INESC-ID, through PEst-OE/EEI/LA0021/2013.

Appendix A. Nominal parameters Acknowledgements See Tables A.1–A.4 The research leading to these results received funding from Fundac¸ão para a Ciência e a Tecnologia via the project Galeno PTDC/SAU-BEB/103667/2008, and from the European Research Council via the Advanced Grant 247035. This work was supported in part by FEDER funds through COMPETE–Programa Operacional Factores de Competitividade and by Portuguese funds through the Center for Research and Development in Mathematics and

Appendix B. Singular values and singular vectors. See Tables B.5–B.8.

34

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

Table B.6 – Singular values and corresponding singular vectors (10−3 ) for results in Fig. 7. SV index SV k1 k2 ˛  C50

1 863

2 342

3.84 −3.04 96.4 −622 777

−256 −203 −934 34.1 144

3 87.2

4 29.6

154 129 −188 −762 −587

5 2.69 − 820 560 93.7 −55.5 −49.8

488 793 −274 169 170

Table B.7 – Singular values (SV) and corresponding singular vectors for results in Fig. 11. The values in the row marked with the superscript * have dimension 10−3 . Index SV

1 3050

2 618

3 383

4 178

5 75.7

6 34.4

7 24.9

8 21.8

9 16.5

10 5.29

k10p k12p k13p k21p k31p ke0p ∗ V1p

−0.7977 −0.0991 −0.3035 −0.1655 0.0109 −0.1232 −0.0126

−0.0398 0.0950 0.0056 −0.0680 0.0046 0.0049 −0.0443

0.3404 −0.4446 0.0076 0.3950 −0.0279 −0.0417 −0.0650

−0.0467 −0.1051 −0.0519 0.0711 −0.0058 −0.0300 −0.0073

0.0676 −0.1479 0.1115 −0.1320 0.0142 −0.0439 −0.0758

0.0294 0.6017 −0.2384 0.1124 −0.0412 0.0415 0.2128

−0.2554 −0.5502 −0.0896 0.2473 0.0009 −0.0127 0.2580

0.1576 0.0159 0.0956 −0.0810 0.0050 −0.2212 0.1923

−0.0881 −0.0987 0.6264 −0.5371 0.1133 0.0742 0.0037

−0.2109 0.0764 0.4993 0.2655 0.0891 −0.0652 0.0029

k10r k12r k13r k21r k31r ke0r V1r C50p C50r  ˇ y0

−0.1883 0.0076 −0.0014 −0.1557 0.0020 −0.0202 −0.0362 −0.3038 −0.0002 0.1972 0.0013 −0.1623

0.0887 −0.0113 0.0030 0.0630 −0.0002 −0.0019 0.9708 −0.0740 −0.0000 −0.0743 0.0004 −0.1388

−0.4231 0.0274 −0.0147 −0.3357 0.0018 −0.0488 0.1210 −0.1641 −0.0004 0.3426 0.0043 −0.2718

−0.1100 0.0097 −0.0028 −0.1168 0.0012 −0.0122 0.1763 −0.0672 −0.0001 0.2172 0.0007 0.9324

0.2211 −0.1777 0.0019 0.6178 −0.0069 −0.0030 −0.0216 −0.1368 −0.0000 0.6761 0.0016 −0.0563

0.0946 −0.0264 0.0146 −0.3357 0.0171 0.0299 0.0228 0.3842 0.0003 0.5343 −0.0063 −0.0827

0.4776 0.0922 0.0075 −0.0666 −0.0040 0.0818 0.0628 0.5569 0.0008 0.0092 −0.0087 −0.0248

0.6098 0.3368 0.0241 −0.3815 0.0031 0.1250 −0.0606 −0.5036 0.0010 0.0756 −0.0125 −0.0060

−0.1935 0.2490 −0.0232 −0.2458 −0.0311 0.0131 0.0339 0.2884 0.0003 0.1922 0.0005 −0.0231

0.1601 −0.7203 −0.0104 −0.2446 −0.0157 −0.0399 −0.0175 −0.0989 −0.0006 −0.0376 0.0063 0.0042

Index SV

11 1.96

12 0.832

13 0.130

14 0.0608

15 0.0313

16 0.0130

k10p k12p k13p k21p k31p ke0p ∗ V1p

0.3028 −0.2400 −0.3782 −0.5847 −0.0659 −0.1901 0.0001

0.0005 0.1086 0.0981 0.0863 0.0044 −0.9377 0.0003

0.0051 −0.0021 0.0019 −0.0160 0.0161 −0.0527 −0.0056

0.0214 0.0223 −0.0665 0.0171 0.2007 −0.0365 −0.0029

0.0121 −0.0020 −0.0288 0.0077 0.2134 0.0123 −0.0234

0.0400 0.0006 −0.0886 0.0149 0.5677 −0.0030 −0.0095

k10r k12r k13r k21r k31r ke0r V1r C50p C50r  ˇ y0

0.0543 −0.4980 0.0031 −0.2478 −0.0019 −0.0052 0.0083 0.1033 −0.0006 −0.0464 0.0061 0.0023

−0.1550 0.0887 0.0149 0.1362 0.0188 0.0427 0.0171 0.1909 0.0011 −0.0474 0.0027 0.0010

0.1236 0.0843 0.1451 −0.0312 0.1209 −0.9664 −0.0089 0.0254 −0.0045 0.0033 0.0536 −0.0001

0.0410 0.0289 − 0.6343 −0.0023 − 0.6950 − 0.1581 −0.0021 0.0023 −0.0107 0.0007 0.1991 0.0000

−0.0145 −0.0008 0.4439 0.0000 − 0.1196 0.1000 0.0002 −0.0001 0.0039 0.0014 0.8551 −0.0000

−0.0222 −0.0057 0.5052 0.0042 − 0.4248 0.0022 0.0006 −0.0000 0.1104 −0.0011 − 0.4685 0.0000

17 0.00343

18 0.00117

19 0.000687

−0.0489 0.0018 0.1058 −0.0071 − 0.7488 0.0005 −0.0796

−0.0001 −0.0001 0.0002 −0.0010 0.0002 0.0007 0.4221

0.0034 −0.0002 −0.0074 −0.0001 0.0539 0.0001 −0.9027

−0.0051 0.0006 0.3319 −0.0026 − 0.5408 −0.0339 −0.0003 0.0003 0.1116 0.0001 −0.0556 −0.0000

0.0028 0.0002 −0.0755 −0.0004 0.0711 −0.0055 −0.0000 −0.0004 0.8988 0.0000 0.0573 0.0000

0.0017 −0.0004 −0.0802 0.0000 0.0901 0.0043 0.0000 −0.0001 0.4092 0.0000 0.0135 0.0000

35

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

Table B.8 – Singular values and corresponding singular vectors (10−3 ) for results in Fig. 14. Index SV d1 d2 l1 l2   C50p C50r m  y0

1 3497 10.9 97.7 10.6 21.2 11.1 52.2 104 −9.08 33.6 −28.8 987

2 724 6.80 60.6 −26.2 −33.6 61.5 357 −46.6 −85.9 117 917 2.22

3 157

4 20.7

5 12.1

6 3.64

24.0 205 −23.2 −24.6 716 408 −177 −104 35 −290 −53.3

15.8 161 −92.3 −115 −640 568 −32.5 −344 185 −253 −49.0

5.59 −9.79 21.7 −2.06 167 526 401 206 −694 −86.2 −48.2

64.5 522 −199 −277 144 −312 151 −605 −310 51.7 −38.7

Index SV

7 1.10

8 0.405

9 0.224

10 0.0202

11 0.00550

d1 d2 l1 l2   C50p C50r m  y0

−87.6 −715 −119 −170 134 71.0 −334 −480 −247 −29.0 110

−52.0 −346 −117 −108 76.3 −48.6 806 −175 401 3.41 −59.8

3.62 7.68 466 759 21.9 6.80 97.7 −439 −49.6 6.70 −35.1

− 163 26.2 − 828 533 1.85 4.01 − 24.9 40.2 − 28.1 −0.232 0.375

− 979 121 141 − 85.7 1.97 0.202 −0.24 −3.58 −5.48 0.267 −0.750

references

[1] American Society of Anesthesiologists, Quality

[2]

[3]

[4]

[5]

[6]

[7]

[8]

Management and Departmental Administration, Continuum of Depth of Sedation: Definition of General Anesthesia and Levels of Sedation/Analgesia, 2009, October. H. Magalhães, J. Lemos, T. Mendonc¸a, P. Rocha, S. Esteves, J. Gaivão, Observer design in switching control of neuromuscular blockade: clinical cases, in: Proc. 28th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC’06), 2006, pp. 5436–5439. G.N. Schmidt, P. Bischoff, T. Standl, G. Lankenau, M. Hilbert, J. Schulte am Esch, Comparative evaluation of narcotrend bispectral index and classical electroencephalographic variables during induction maintenance and emergence of a propofol/remifentanil anesthesia, Anesthesia & Analgesia 98 (5) (2004) 1346–1353. J.W. Johansen, P.S. Sebel, Development and clinical application of electroencephalographic bispectrum monitoring, Anesthesiology 93 (2000) 1336–1344. A. Raue, C. Kreutz, T. Maiwald, U. Klingmüandller, J. Timmer, Addressing parameter identifiability by model-based experimentation, IET Systems Biology 5 (2) (2011) 120–130. S. Audoly, L. D’Angio, M. Saccomani, C. Cobelli, Global identifiability of linear compartmental models-a computer algebra algorithm, IEEE Transactions on Biomedical Engineering 45 (1) (1998) 36–47. S. Audoly, G. Bellu, L. D’Angio, M. Saccomani, C. Cobelli, Global identifiability of nonlinear models of biological systems, IEEE Transactions on Biomedical Engineering 48 (1) (2001) 55–65. L.Y. Wang, G. Yin, H. Wang, Reliable nonlinear identification in medical applications, in: Proc. 13th IFAC Symposium on System Identification (SYSID’03), 2003.

[9] I. Gustavsson, L. Ljung, T. Söderström, Survey paper: identification of processes in closed loop – identifiability and accuracy aspects, Automatica 13 (1977) 59–75. [10] T. Söderström, P. Stoica, System Identification, Prentice-Hall, Hemel Hempstead, UK, 1989. [11] J.A. Jacquez, T. Perry, Parameter estimation: local identifiability of parameters, American Journal of Physiology – Endocrinology And Metabolism 258 (4) (1990) E727–E736. [12] K. Soltesz, K. van Heusden, G.A. Dumont, T. Hägglund, C.L. Petersen, N. West, J.M. Ansermino, Closed-loop anesthesia in children using a pid controller: a pilot study, in: Proc. IFAC Conference on Advances in PID Control (PID’03), 2012. [13] J. Lemos, J. Gomes, B. Costa, T. Mendonc¸a, A. Coito, Batch identification of neuromuscular blockade models, in: Proc. 19th Mediterranean Conference on Control Automation (MED’11), 2011, pp. 503–508. [14] J. Lemos, J. Gomes, B. Costa, P. Gambus, E. Jensen, T. Mendonc¸a, A. Coito, A bi-phase algorithm to identify hypnotic models of patients subject to deep sedation for ultrasonographic endoscopy, in: Proc. 16th IFAC Symposium on System Identification (SYSID’12), 2012, pp. 775–780. [15] M.M. Silva, T. Wigren, T. Mendonc¸a, Nonlinear identification of a minimal neuromuscular blockade model in anesthesia, IEEE Transactions on Control Systems Technology 20 (1) (2012) 181–188. [16] M.M. Silva, T. Mendonc¸a, T. Wigren, Online nonlinear identification of the effect of drugs in anaesthesia using a minimal parameterization and BIS measurements, in: Proc. American Control Conference (ACC’10), 2010, pp. 4379–4384. [17] T. Wigren, Recursive prediction error identification method using the nonlinear Wiener model, Automatica 29 (4) (1993) 1011–1025. [18] B. Weatherley, S. Williams, E. Neill, Pharmacokinetics, pharmacodynamics and dose–response relationships of atracurium administered i.v., British Journal of Anaesthesia 55 (1983) 39s–45s.

36

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 23–36

[19] T. Mendonc¸a, P. Lago, PID control strategies for the automatic control of neuromuscular blockade, Control Engineering Practice 6 (10) (1998) 1225–1231. [20] K. Godfrey, Compartmental Models and their Application, Academic Press, Oxford/Great Britan, 1983. [21] M.M. Silva, Prediction error identification of minimally parameterized wiener models in anesthesia, in: Proc. 18th IFAC World Congress, 2011, pp. 5615–5620. [22] M.M. Silva, R. Rabic¸o, T. Mendonc¸a, T. Wigren, Control of rocuronium-induced neuromuscular blockade via online identification of a two-parameters wiener model, in: Proc. 16th IFAC Symposium on System Identification (SYSID’12), 2012, pp. 571–576. [23] C.F. Minto, T.W. Schnider, T.G. Short, K.M. Gregg, A. Gentilini, S.L. Shafer, Response surface model for anesthetic drug interactions, Anesthesiology 92 (2000) 1603–1616. [24] E. Gepts, F. Camu, I. Cockshott, E. Douglas, Disposition of propofol administered as constant rate intravenous infusions in humans, Anesthesia & Analgesia 66 (12) (1987 December) 1256–1263. [25] C. Minto, T. Schnider, T. Egan, E. Youngs, H. Lemmens, P. Gambus, V. Billard, J. Hoke, K. Moore, D. Hermann, K. Muir, J. Mandema, S. Shafer, Influence of age and gender on the pharmacokinetics and pharmacodynamics of remifentanil. I. Model development, Anesthesiology 86 (1997) 10–23. [26] A. Gentilini, C. Frei, A. Glattfedler, M. Morari, T. Sieber, R. Wymann, T. Schnider, A. Zbinden, Multitasked closed-loop control in anesthesia, IEEE Engineering in Medicine and Biology Magazine 20 (1) (2001) 39–53. [27] M.M. Silva, T. Wigren, T. Mendonc¸a, Nonlinear identification of the human response to anaesthesia by the use of BIS

[28]

[29]

[30]

[31]

[32]

[33]

measurements, Tech. Rep. 2010-008, Dept. Information Technology, Uppsala University, Uppsala, Sweden, 2010, March), Available: http://www.it.uu.se/research/publications/reports/ S. Vinga, K. Thomaseth, J. Lemos, A. Neves, H. Santos, A. Freitas, Structural analysis of metabolic networks: a case study on lactococcus lactis, in: Proc. 8th Portuguese Conference on Automatic Control (CONTROLO’2008), 2008, pp. 4379–4384. P. Lago, T. Mendonc¸a, L. Gonc¸alves, On-line autocalibration of a PID controller of neuromuscular blockade, in: Proc. IEEE International Conference on Control Applications, 1998, pp. 363–367. T.W. Schnider, C.F. Minto, P.L. Gambus, C. Andresen, D.B. Goodale, S.L. Shafer, E.J. Youngs, The influence of method of administration and covariates on the pharmacokinetics of propofol in adult volunteers, Anesthesiology 8 (1998) 1170–1182. A.R. Absalom, M.M.R.F. Struys, Overview of Target Controlled Infusions and Total Intravenous Anaesthesia, 2nd ed., Academia Press, Gent, 2007. M. Mertens, E. Olofsen, F. Engbers, A. Burm, J. Bovill, J. Vuyk, Propofol reduces perioperative remifentanil requirements in a synergistic manner: response surface modeling of perioperative remifentanil–propofol interactions, Anesthesiology 99 (2) (2003) 347–359. B. Marsh, M. White, N. Morton, G.N.C. Kenny, Pharmacokinetic model driven infusion of propofol in children, British Journal of Anaesthesia 67 (1) (1991) 41–48.