Latin American Women Filmmakers: Production, Politics, Poetics

Pharmacodynamic modeling of propofolinduced general anesthesia in young adults The MIT Faculty has made this article openly available. Please share how this access benefits you. Your story matters. Citation Chakravarty, Sourish et al. "Pharmacodynamic modeling of propofol-induced general anesthesia in young adults." 2017 IEEE Healthcare Innovations and Point of Care Technologies (HI-POCT), November 2017, Bethesda, Maryland, USA, Institute of Electrical and Electronics Engineers (IEEE), December 2017. © 2017 IEEE As Published https://dx.doi.org/10.1109/hic.2017.8227580 Publisher Institute of Electrical and Electronics Engineers (IEEE) Version Author's final manuscript Citable link https://hdl.handle.net/1721.1/123300 Terms of Use Creative Commons Attribution-Noncommercial-Share Alike Detailed Terms https://creativecommons.org/licenses/by-nc-sa/4.0/ Pharmacodynamic modeling of propofol-induced general anesthesia in young adults Sourish Chakravarty ∗ , Ksenia Nikolaeva ‡ , Devika Kishnan ‡ , Francisco J. Flores Patrick L. Purdon † ¶ , Emery N. Brown ∗ † ‡ § ¶ ∗ ¶ k, ∗ The Picower Inst. for Learning and Memory, MIT, Cambridge, MA; † Brain and Cognitive Sci., MIT, Cambridge, MA; ‡ Inst. of Med. Engin. and Sci. , MIT, Cambridge, MA; § MIT-Harvard Hlth. Sci. and Technol., Cambridge, MA; ¶ Anesthesia, Critical Care and Pain Medicine, Mass. Gen. Hosp., Boston, MA; k Harvard Medical School, Boston, MA; Abstract—Target controlled infusion (TCI) of intraveneous anesthetics can assist clinical practitioners to provide improved care for General Anesthesia (GA). Pharmacokinetic/Pharmacodynamic (PK/PD) models help in relating the anesthetic drug infusion to observed brain activity inferred from electroencephalogram (EEG) signals. The parameters in popular population PK/PD models for propofol-induced GA (Marsh and Schnider models) are either verified based on proprietary functions of the EEG signal which are difficult to correlate with the neurophysiological models of anesthesia, or the marker itself needs to be estimated simultaneously with the PD model. Both these factors make these existing paradigms challenging to apply in real-time context where a patient-specific tuning of parameters is desired. In this work, we propose a simpler EEG marker from frequency domain description of EEG and develop two corresponding PK/PD modeling approaches which differ in whether they use existing population-level PK models (approach 1) or not (approach 2). We use a simple deterministic parameter estimation approach to identify the unknown PK/PD model parameters from an existing human EEG data-set. We infer that both approaches 1 and 2 yield similar and reasonably good fits to the marker data. This work can be useful in developing patient-specific TCI strategies to induce GA. Index Terms—General Anesthesia, Target Controlled Infusion I. I NTRODUCTION General Anesthesia (GA) is a vital element of surgical procedures. Clinical practitioners of anesthesia introduce hypnotic drugs into a patient during a major surgery with the objective to achieve a state of GA - a reversible state of unconciousness where the patient does not move or feel pain or have memory of the surgical procedure while physiological stability is maintained. [1]. Monitoring the electroencephalogram (EEG) [2] can provide useful additional information to the anesthesiologist for inferring the brain’s state under GA. In this context, an integrated pharmacokinetic (PK) and EEGbased pharmacodynamic (PD) modeling framework provides a principled approach to model, respectively, changes in the amount of the drug in the therapeutic site and how these changes affect an EEG-derived marker indicating the level of GA. This PK/PD paradigm of describing the dose-effect relationship has found use in Target Controlled Infusion (TCI) systems that are designed to change the dose automatically based upon prescribed target plasma or effect site concentrations of the drug (anesthetic agent) [3]. For propofol-induced GA, Marsh [4] and Schnider [5] models are two common PK models used in TCI systems to describe the drug concentration dynamics in the blood/ plasma [6]. Both PK models have also been used to link the estimated drug concentration dynamics in the brain to that of a prescribed EEG-marker, Bispectral Index (BIS) [7] for Marsh and Canonical Univariate Parameter (CUP) [8] for Schnider, via respective PD models [6]. BIS, a depth of anesthesia (DOA) marker, is an empirically established number combining both frequency and time domain features of the EEG signal [9] and, being proprietary, its operational definition is not publicly available. The CUP marker is a weighted sum of logarithms of band-wise power in the instantaneous power spectral density (PSD). The mechanism underlying such band-wise dynamics in the PSD under anesthesia can potentially be reasoned using network models of neurons. For example, appearance of spatially coherent frontal alpha (8 12 Hz) oscillations in the EEG [10] observed simultaneous with loss of conciousness (LOC) during propofol anesthesia could be explained by existing neurophysiological model of thalamo-cortical circuits [11]. However, the Schnider PK/PD model was developed by simultaneously estimating both the weights of the CUP marker and the parameters of the corresponding PD model [8]. The nature of coupling between the definitions of the marker and the functional form of the PD model makes the verification of the Schnider PK/PD model challenging since refinement of the PK/PD model would go hand-in-hand with re-estimation of CUP marker weights. Thus, there is a need for spectral markers that can be defined independent of a PK/PD model, and are able to indicate the level of of GA based on known features of EEG oscillations characteristic of propofol-induced GA [10]. In this work we first simplify the Schnider CUP marker by redefining a spectral marker using powers in fixed frequencybands of the PSD which reflect the known signatures associated with propofol-induced GA [10]. We model this marker using a bi-phasic PD model similar to one used by Schnider et. al. [8].We propose and demonstrate two PK/PD modeling approaches and corresponding parameter estimation schemes that can apply our choice of marker to human EEG data reflective of anesthetic induction with propofol. II. M ETHODS A. Data collection This study was approved by the Human Research Committee at Massachusetts General Hospital and involved 10 healthy young adult volunteers aged between 18-36 years undergoing induction into and recovery from GA [10]. For each subject, intravenous infusion of propofol was administered so as to maintain a target effect site concentration at five monotonically increasing levels where each level was maintained for 14 minutes during the induction phase (based on Schnider PK/PD model [5], [8]). A similar step-wise infusion was planned to effect a gradual emergence. In this paper we focus on the data collected in the induction phase only. The time points for LOC and return of conciousness were determined based on each subject’s response to auditory stimuli. Subjects were instructed to keep their eyes closed during the study sessions. EEG was recorded continuously from each subject using a 64-channel BrainVision MRI Plus System (Brain Products) with a sampling rate of 5 kHz, resolution 0.5 µV . As part of the pre-processing, EEG signals were passed through an anti-aliasing filter and downsampled to 250 Hz. A nearest-neighbor Laplacian referencing approach was used to minimize the spatial blurring. The current work uses such pre-processed EEG data from one of the frontal electrodes. Fig. 1. (a, b, c) Infusion rates of propofol, (d, e, f) Multitaper Spectrogram estimate with color scale in dB, (g, h, i) Marker observation (blue dots) and smoothed estimated observation (blue solid line) for 3 subjects (A, B, C). The relevant parameters chosen for MT analysis are: time-window of 10 s with overlap of ∆ = 5 s between consecutive windows, time-halfbandwidth product of 5 and 4 tapers. C. Pharmacodynamic (PD) model of the chosen marker To model the bi-phasic trend of the chosen EEG marker in Eq. (1) with monotonic changes in anesthetic level, we choose a pharmacodynamic model with two sigmoid functions connected at a junction point [8] described as (1) (1) R1γ1 / (1 + R1γ1 ) , if 0 < x(t) ≤ xh , yt = E0 − Emax (2) B. EEG Marker from spectrogram For each subject, a multitaper (MT) spectrogram [12] is obtained from the EEG time-series using Chronux toolbox [13] in MATLAB 2017a software package. We calculate powers in slow (0.1 - 1 Hz), delta (1 - 4 Hz) and beta (14 30 Hz) bands from MT PSD estimates, respectively denoted (s) (d) (b) as Pt , Pt and Pt , at time instant t. The chosen marker is, yt⋆ = (s) 10 log10 ((Pt + (d) (b) Pt )/Pt ) . (1) For propofol induction in the 10 subjects the overall trend in the marker yt⋆ displays an initial reduction (relative to a baseline value) for a brief duration followed by a gradual rise which asymptotically approaches an upper bound. The initial dip in the marker under low doses of propofol is due to increased beta activity [10] and can potentially be explained by network models of neurons [14]. The subsequent rise after LOC is due to reduction in beta activity and rise in low frequency power [10]. Prior to model fitting, the marker observation to-be-fitted is estimated by applying a randomwalk fixed interval smoother. The relevant parameters for the filtering and smoothing operations are prescribed as follows: initial value and variance of scalar hidden state is assumed to be same as the mean and standard error of mean from a 10 minute segment of {yt⋆ } prior to start of drug infusion. The variances of the observation and state (zeromean Gaussian white) noises are respectively taken as 1× and 0.001× variance of this initial data segment. For the ensuing discussion we shall indicate the smoothed estimate also as yt⋆ . Fig. 1 illustrates the quality of the smoothing step for 3 subjects. (2) or, = E0 + Emax R2γ2 / (1 + R2γ2 ) , if xh < x(t) , (1) (2) (2) where, R1 = x(t)/x50 and R2 = (x(t) − xh )/(x50 − xh ) (1) (2) (1) (2) and Emax , Emax , γ1 , γ2 , x50 , x50 and xh are positive realvalued scalars. The term x(t) can be regarded as a positive continuous function of the infusion history that indicates the drug-induced level of anesthesia. A higher value of x(t) would indicate a deeper level of GA. Additionally, we assume (2) E0 (1) (1) (1) (1) (xh /x50 )γ1 /(1 + (xh /x50 )γ1 ) . (3) = E0 − Emax The marker observation yt⋆ is modeled using yt which shows (2) (1) a decrease from a baseline value of E0 to E0 with rise in x(t) upto a critical cutoff level xh and for values of x(t) > xh this marker rises and then approaches a maximum value of (2) (2) E0 + Emax asymptotically for x(t) ≫ xh . Two approaches are proposed to model x(t) for our chosen marker. 1) Approach 1: In this approach, we assume that a population PK model is available to describe the time-course of drug amount in the blood/plasma (denoted by xc (t)) for a known drug infusion rate {u(t)}. In this work we use the three-compartment Schnider PK model [15] as the known population PK model to pre-calculate xc (t) for a subject. To model the drug dynamics in the brain, an additional effectsite compartment of very small size is conceived to represent the brain [8]. This effect-site compartment is assumed to exchange drug with the plasma based on first-order diffusion kinetics. We interpret x(t) from Eq. (2) as amount of drug in the effect-site (denoted by xe (t)) which can be calculated from xc (t) as, Zt (4) x(t) ≡ xe (t) = kce e−ke0 (t−τ ) xc (τ )dτ . 0 Based on the assumption that addition of this effect compartment does not affect the dynamics of the original xc (t), one can arbitrarily set kce = ke0 /10000 [15]. Thus only one parameter in Eq (4) needs to be estimated in addition to the PD model parameters in Eq. (2). Our approach 1, is in principle similar to Schnider et. al. [8], but this earlier work also estimated the weights in the CUP marker along with the PD model parameters and ke0 . The approach 1 serves two purposes in our work: (1) It enables us to submit the 10 patient data-set considered here to existing tools available in PK/PD modeling literature, and (2) helps in benchmarking the quality of the modeling using the newer approach 2. Approach 1 is generalizable to the case when the underlying known Schnider PK model (for xc (t)) is replaced by the Marsh model. 2) Approach 2: In our second approach, we consider a more general scenario where no population PK model is available. In this approach we use the mathematical representation of a two-compartmental PK (PK2) model to describe x(t) as a non-dimensional scalar that can also be interpreted as an index of the level of anesthesia. A PK2 model is described by a bivariate state vector that evolves as,        1 x1 −(k12 + k10 ) k21 ẋ1 + u(t) (5) = k12 −k21 x2 ẋ2 0 bounds. For both approaches, the parameter estimation problem is stated as the following optimization problem yielding an optimum value of parameter set Θ as, e = argmin F (Θ) ; F (Θ) = F1 (Θ) + F2 (Θ) + F3 (Θ) Θ Θ F1 = F2 = F3 = X tk ∈τ (2) (E0 (2) (E0 yt k −
2 yt⋆k (7) ⋆ /(2Nt (ymax 2 ⋆ ymin ) ) − ⋆ ⋆ ⋆ − ymin )2 /2(ymax − ymin )2 + (2) Emax − ⋆ ⋆ ymax )2 /2(ymax − (8) (9) ⋆ ymin )2 (10) where, τ = {tk = k∆, ∀ k ∈ [0, Nt − 1]} and Nt is the length of smoothed observation sequence {yt⋆ } and t0 = 0 denotes the time-point when infusion was started. For approach 1, the objective function is minimized over With intial conditions given as x1 (0) = 0 and x2 (0) = 0, the analytic solution for x2 (t) is given as, k12 x2 (t) = λ1 − λ 2 Zt  0  e(t−τ )λ1 − e(t−τ )λ2 u(τ ) dτ (6) where, λ1 and λ2 are assumed to be real and distinct eigenvalues of the rate-constant matrix such that λ2 < λ1 < 0. If one uses x2 (t) as the x(t) in Eq. (2), then the multiplying factor outside the integral in Eq. (6) will be multiplied with (1) (2) x50 , x50 and xh in Eq. (2). Therefore, one can simplify the problem further by considering x(t) ≡ x2 d/k12 where, d ≡ λ1 − λ2 > 0. The unkown PK parameters are λ1 and d. D. Parameter estimation In the PK/PD formulation based on Eqs.(2), (4) and (6) the unknown model parameter sets are Θ(1) = (1) (2) {ke0 , Emax , Emax , γ1 , γ2 , α1 , α2 , αh } for approach 1 and (1) (2) Θ(2) = {λ1 , d, Emax , Emax , γ1 , γ2 , α1 , α2 , αh } for ap(1) proach 2. We calculate E0 by taking mean over an initial segment of baseline smoothed marker data (before drug infusion). To aid the optimization program, we identify the ⋆ ⋆ minimum (ymin ) and maximum (ymax ) values of marker data during the induction phase. It is assumed here the trajectory u(t) is such that it yields a x(t) which has a monotonically non-decreasing trend. This assumption is compatible to the data set considered here [10] and the sequence x(t) has a global maximum, xmax , in the domain of interest. The (1) following re-parametrization: xh = αh xmax , x50 = α1 xmax , (2) x50 = α2 xmax yields non-dimensional shape parameters αh , α1 and α2 which can then be constrained within reasonable Fig. 2. (a, b, c) Model fits from two approaches overlaid on the observed marker, (d, e, f) Normalized x(t) with junction points of two sigmoids indicated for 3 subjects (A, B, C). The pair of optimal values e (1) ), F (Θ e (2) )) are: (0.1168, 0.1411) × 10−3 for A, (0.5137, (F (Θ 0.6331) × 10−3 for B, (1.4472, 1.3641) × 10−3 for C. a uniformly spaced grid ke0 ∈ [0.01, 0.2] min−1 such that for each choice of ke0 on this search grid, a local gradientbased optimization problem is solved for the rest of the parameters Θ(1) . For a given choice of ke0 and {u(t)}, we set initial guess αh(o) as x(th )/xmax where th is the time-point ⋆ corresponding to yt⋆h = ymin . Here, we use the convention that θ(o) will denote the initial guess for a parameter θ. We found our choice of domain of ke0 to be suitable for the current data-set based on intial parameter sweeps over wider ranges of ke0 that included ke0 = 0.456 min−1 (as reported in Schnider et. al. [8]). For approach 2, analogous to the optimization procedure in approach 1, the objective function is minimized over a uniformly spaced grid λ1 ∈ [−0.2, −0.01] min−1 such that for each choice of λ1 in this search grid a local gradient-based optimization is solved for the rest of the parameters in Θ(2) . For a given choice of λ1 , we set d(o) = −λ1 , αh(o) = 0.3 and use 0.001 ≤ d ≤ 3 as an additional constraint. Common to both approaches, the rest of the optimization problem prescription is as follows: (1) (1) ⋆ Emax(o) = 2(E0 − ymin ), γ1(o) = 1, α1(o) = 0.1, (2) ⋆ ⋆ Emax(o) = ymax − ymin , γ2(o) = 1, α2(o) = 0.5 + αh(o) /2, (1) (1) (1) 0.8Emax(o) ≤ Emax ≤ 100Emax(o) , 0.5 ≤ γ (1) ≤ 6, 0.001 ≤ α1 ≤ 2, 1 ≤ γ2 ≤ 6, 0.8αh(o) ≤ α2 ≤ 1, 0.1 ≤ (2) (2) (2) αh ≤ 0.8, 0.8Emax(o) ≤ Emax ≤ 1.2Emax(o) and αh < α2 . Both optimization problems are solved in MATLAB 2017a using Sequential Quadratic Programming. III. R ESULTS We apply both the approaches, 1 & 2, to model the estimated marker from all 10 subjects. We identify 3 subjects, A, B and C, respectively corresponding to minimum (best fit), median and maximum (worst fit) among 10 values of e (1) ) (also the same subjects in Fig. 1). The fits for these F (Θ 3 subjects from both approaches are presented in Fig. 2(a, b, c). Both approaches yield near-identical fits to the markers from the current data-set. Although, the exact values of x(t) are different from two approaches, the shape of the curves from both approaches are qualitatively similar Fig. 2(d, e, f) to each other. The spread of the parameters estimated well on the current data-set. Thus further research along the lines of Approach 2 might be fruitful. Additionally, simplicity of the marker definition and the setup of the optimization problem enables the use of similar parameter estimation schemes for both approaches. This helps in establishing an easy-to-implement process to benchmark future developments based on approach 2. Also, based on this retrospective analysis, we conclude that our choice of initial guess and optimization problem formulation yields a local minimum within the parametric domain considered here. This can potentially help in fast real-time estimation of patient-specific parameters using marker data from induction phase. Acknowledgments: This work was partially supported by National Institutes of Health (NIH) Transformative Research Award R01-GM104948 (to E.N.B.), by NIH Award P01-GM118629 (to E.N.B,), by funds from Massachusetts General Hospital (to E.N.B.), by funds from the Picower Institute for Learning and Memory (to E.N.B.). R EFERENCES Fig. 3. Box-plots from Approach 1 (a) and 2 (b) for all 10 subjects. Each parameter field is normalized by the median value for better visualization. The respective values for the median corresponding to optimum parameters (row-wise from left to right on the xaxes) for subfigure (a) are {0.0005, 0.0350, 7.6747, 13.6887, 0.8556, 2.5106, 2.0000, 0.3691, 0.1660} and for subfigure (b) are {0.0006, −0.0100, 2.2674, 7.6581, 13.6559, 0.8972, 2.5315, .1.0000, 0.4162, 0.1734} from all subjects using both approaches can be visualized via the box-plots in Fig. 3. To compare the estimates of a common parameter θe between Figs. 3(a) and (b), we perform a two-sample t-test and obtain 95% confidence intervals (CI) for the difference between θe(1) − θe(2) . These are given as (−3.3293, 3.3516), (−3.4477, 3.4586), (−0.5686, 0.4320), (−0.6365, 0.4932), (−0.1594, 1.2259), (−0.1338, 0.1023) (1) and (−0.1124, 0.1160) for θ respectively indicating E0 − (2) (2) E0 , Emax , γ1 , γ2 , α1 , α2 and αh . Based on this analysis, we infer that the difference in the means for all the common parameters from the two approaches are not significantly different. IV. D ISCUSSION Our work expands the previous work from Schnider et. al. [5] on marker and PD model for propofol-induced LOC. We simplified the definition of the CUP marker so as to make it independent of the PD modeling without losing the bi-phasic signature characteristic of propofol. Furthermore, we utilized a conventional PK/PD framework of approach 1 to model the marker dynamics and also estimate the underlying drug-induced level of anesthesia x(t). With the connection established between our current data-set and classical PK/PD modeling framework of propofol, we follow the more general approach 2 wherein we view x(t) as a latent state within a two-compartment PK-inspired modeling framework. The latter procedure can be applied even when a population PK model is unavailable. For the same drug infusion and chosen electrode-location, we observe that both approaches are able to fit the marker data reasonably [1] E. N. Brown, R. Lydic, and N. D. Schiff, “General anesthesia, sleep, and coma,” New England Journal of Medicine, vol. 363, no. 27, pp. 2638–2650, 2010. [2] I. J. Rampil, “A primer for eeg signal processing in anesthesia,” Anesthesiology: The Journal of the American Society of Anesthesiologists, vol. 89, no. 4, pp. 980–1002, 1998. [3] A. R. Absalom and K. P. Mason, Total Intravenous Anesthesia and Target Controlled Infusions. Springer, 2007. [4] B. Marsh, M. White, N. Morton, and G. Kenny, “Pharmacokinetic model driven infusion of propofol in children,” BJA: British Journal of Anaesthesia, vol. 67, no. 1, pp. 41–48, 1991. [5] T. W. Schnider, C. F. Minto, P. L. Gambus, C. Andresen, D. B. Goodale, S. L. Shafer, and E. J. Youngs, “The influence of method of administration and covariates on the pharmacokinetics of propofol in adult volunteers,” Anesthesiology: The Journal of the American Society of Anesthesiologists, vol. 88, no. 5, pp. 1170–1182, 1998. [6] A. Absalom, V. Mani, T. De Smet, and M. Struys, “Pharmacokinetic models for propofol - defining and illuminating the devil in the detail,” British journal of anaesthesia, vol. 103, no. 1, pp. 26–37, 2009. [7] M. M. Struys, T. De Smet, B. Depoorter, L. F. Versichelen, E. P. Mortier, F. J. Dumortier, S. L. Shafer, and G. Rolly, “Comparison of plasma compartment versus two methods for effect compartment– controlled target-controlled infusion for propofol,” Anesthesiology, vol. 92, no. 2, pp. 399–406, 2000. [8] T. W. Schnider, C. F. Minto, S. L. Shafer, P. L. Gambus, C. Andresen, D. B. Goodale, and E. J. Youngs, “The influence of age on propofol pharmacodynamics,” The Journal of the American Society of Anesthesiologists, vol. 90, no. 6, pp. 1502–1516, 1999. [9] J. C. Sigl and N. G. Chamoun, “An introduction to bispectral analysis for the electroencephalogram,” Journal of Clinical Monitoring and Computing, vol. 10, no. 6, pp. 392–404, 1994. [10] P. L. Purdon, E. T. Pierce, E. A. Mukamel, M. J. Prerau, J. L. Walsh, K. F. K. Wong, A. F. Salazar-Gomez, P. G. Harrell, A. L. Sampson, A. Cimenser et al., “Electroencephalogram signatures of loss and recovery of consciousness from propofol,” Proceedings of the National Academy of Sciences, vol. 110, no. 12, pp. E1142–E1151, 2013. [11] S. Ching, A. Cimenser, P. L. Purdon, E. N. Brown, and N. J. Kopell, “Thalamocortical model for a propofol-induced α-rhythm associated with loss of consciousness,” Proceedings of the National Academy of Sciences, vol. 107, no. 52, pp. 22 665–22 670, 2010. [12] D. J. Thomson, “Spectrum estimation and harmonic analysis,” Proceedings of the IEEE, vol. 70, no. 9, pp. 1055–1096, 1982. [13] H. Bokil, P. Andrews, J. E. Kulkarni, S. Mehta, and P. P. Mitra, “Chronux: a platform for analyzing neural signals,” Journal of neuroscience methods, vol. 192, no. 1, pp. 146–151, 2010. [14] M. M. McCarthy, E. N. Brown, and N. Kopell, “Potential network mechanisms mediating electroencephalographic beta rhythm changes during propofol-induced paradoxical excitation,” Journal of Neuroscience, vol. 28, no. 50, pp. 13 488–13 504, 2008. [15] S. L. Shafer and K. M. Gregg, “Algorithms to rapidly achieve and maintain stable drug concentrations at the site of drug effect with a computer-controlled infusion pump,” Journal of Pharmacokinetics and Pharmacodynamics, vol. 20, no. 2, pp. 147–169, 1992.