Real-Time Cubature Kalman Filter Parameter Estimation of Blood Pressure Response Characteristics Under Vasoactive Drugs Administration

Real-T ime Cubature Kalman Filter Parameter Estimation of Blood Pressure Response Characteristics Under V asoacti ve Drugs Administration* Shahin T asoujian, Saeed Salav ati, Karolos Grigoriadis, and Matthew Franchek Abstract — Mathematical modeling and real-time dynamics identification of the mean arterial blood pressur e (MAP) response of a patient to vasoactive drug infusion can provide a reliable tool for automated drug administration and there- for e, reduce the emergency costs and significantly benefit the patient’ s MAP regulation in an intensive care unit. T o this end, a dynamic first-order linear parameter -varying (LPV) model with varying parameters and varying input delay is considered to capture the MAP response dynamics. Such a model effectively addresses the complexity and the intra- and inter-patient variability of the physiological response. W e discretize the model and augment the state vector with model parameters as unknown states of the system and a Bayesian- based multiple-model square r oot cubature Kalman filtering (MMSRCKF) approach is utilized to estimate the model time- varying parameters. Since, unlike the other model parameters, the input delay cannot be captur ed by a random-walk process, a multiple-model module with a posterior probability estimation is implemented to provide the delay identification. V alidation results confirm the effectiveness of the proposed identification algorithm both in simulation scenarios and also using animal experiment data. I . I N T RO D U C T I O N Fast-acting vasoacti ve medications are often used as a potentially vital medical intervention to address the patients’ hemodynamics instability and regulate the blood pressure to a desired target v alue in numerous clinical and emergency resuscitation procedures. The vasoacti ve drugs are di vided into two main categories: (1) v asodilator drugs that are being administered to treat individuals with hypertension. This type of drug widens blood vessels, and thereby , can rapidly lower the blood pressure. One of the most effecti ve vasodilator drugs is sodium nitroprusside (SNP) which has been used to treat elev ated blood pressure in various clinical scenarios such as post-surgical care, childbirth, and treating the typical high blood pressure disorders [1], [2]. (2) V aso- pressor drugs such as phenylephrine (PHP), v asopressin, and norepinephrine increase the blood pressure by stimulating the depressed cardiovascular system causing vasoconstric- tion. These medicines are being used to treat patients with hypotensi ve symptoms in different medical scenarios such as hemorrhage, traumatic brain injury , and septic shock [3]. Recently , the automated closed-loop administration of va- soactiv e drugs for the mean arterial pressure (MAP) control and regulation purposes, has gained significant attention * This paper is a preprint of a paper submitted to the IEEE American Control Conference 2020. S. T asoujian, S. Salavati, K. Grigoriadis, and M. Franchek are with the Department of Mechanical Engineering at the University of Houston, Houston, TX 77004 USA (e-mail: stasoujian, ssalavatidezfuli, karolos, mfranchek@uh.edu). in clinical care [2], [4]. The automated dosage adminis- tration procedure utilizes the wealth of feedback control and surpasses the manual drug infusion with a syringe, in terms of accuracy , timeliness, cost ef ficiency , and reliability . Nev ertheless, in order to have an accurate automated drug administration and to be able to implement robust model- based control strategies, an explicit dynamical model that accurately describes MAP response behavior to the drug infusion is needed. Howe ver , based on experiments, there are significant variations in the patient’ s physiological response to the drug infusion [5]. This pharmacological variation of the patient’ s MAP response to the drug causes the model parameters to vary ov er time for an individual, as well as, from patient-to-patient. Consequently , due to such intra- and inter-patient physiological variations, a mathematical model with fixed parameters cannot be adequate to capture the patient’ s complex MAP dynamics. Moreover , although a robust controller design can guarantee the stability and target MAP tracking of the patient’ s v arying closed-loop system, other properties such as the settling time, rising time, damping characteristics, and disturbance rejection can significantly degrade when the patients parameters are off the control parameters. In this regard, a real-time MAP response modeling and efficient parameter estimation scheme is essential. There are multiple approaches for the estimation of dy- namical system parameters. In the first approach, a lin- earization of the dynamical system is used to perform the parameter estimation. Extended Kalman filtering (EKF) is one of the widely used methods. Howe ver , it is only ap- plicable to systems with mild nonlinearities, and it requires the Jacobian matrix computation. Moreover , numerical errors due to truncation and con vergence problems are likely in EKF and other local approximation based estimators [6]. In the other approach, known as the sampling method, the nonlinear representation of a system is used to estimate the parameters via a filter such as the unscented Kalman filter (UKF) which leads to more accurate estimation. In UKF , a set of weighted sampling points propagates via the nonlinear function of the system. Howe ver , for higher-order systems, UKF is prone to numerical instability since the weights of the sigma points may become negati ve [7]. Another sampling method is particle filtering (PF) which is an iterativ e Monte Carlo based method to compute the posterior probability distribution of the state of a nonlinear system even with non- Gaussian noise. PF requires a large set of randomly generated particles to approximate the posterior probability density function. Under an increase in the number of iterations, PF encounters particle degradation and depletion. In order to ov ercome such issues, the authors in [6] have proposed a Bayesian filtering framework known as cubature Kalman filtering (CKF). The sample points in the CKF algorithm propagate via equally valued cubature points which are twice the size of the system nonlinear function. It uses a spherical- radial cubature rule to generate the weighted sum of sampling points to approximate the integrals in Bayesian estimation. CKF demonstrates better nonlinear performance, stability , and accuracy compared to EKF , PF , and UKF [7]. In [8], the authors have used a first-order model with delayed measurements to describe the MAP dynamics in response to hypotensiv e drugs. They have pre-identified the parameters using dose-response characteristics induced by a rectangular test signal while trying to av oid any adverse effect on the patient. If the identified parameters are not within the prescribed bounds, then the experiment will be repeated. Nonetheless, in the case of outlying identification results, the worst-case parameters are used. The output is filtered by a constant filter which has been derived using trial and error . The delay is determined via the response settling characteristics. In another work, a generalized fuzzy neural network framework has been studied for the estimation and control of MAP dynamics in response to vasodilator drugs [9]. The parameters have been assumed to be nonlinear functions of the measured MAP . This method requires a training dataset and an ef fective learning algorithm for the ar - tificial scheme. Moreov er , ov erparameterization and a proper number of perceptrons remain as other obstacles. Based on the time delay model introduced in [10], the authors in [11] have proposed discrete-time parameters update laws. Howe ver , the procedure of the parameters identification of the original model has not been addressed. A bank of Kalman filters (KFs) augmented with a posterior probability estimator to match a candidate model to that of the patient has been designed in [12]. Each KF is responsible for generating the state vector updates for the next step, and the Kalman gain is assumed to be generated a priori. Then, by calculating the residual of the actual and generated output, the state vector is updated accordingly . In order to capture the varying time delay , which a conv entional KF is not capable of, the multiple-model (MM) approach has been adopted through which fi ve equally spaced delay blocks from 10 s to 50 s , each is considered to be cascaded with the same bank of KFs. The recursiv e posterior probability estimation is calculated for each residual to update the input with the most likely delay . In a similar approach, [2] has examined KF for the estimation of the MAP dynamics parameters in hypertension. In this work, authors have utilized the model introduced first in [10]. They hav e discretized and transformed the infinite- dimensional model into a linear one that accommodates an input with three backward steps. Then, the parameters are gathered in a vector that is updated through the KF approach. The control parameters are updated using heuristic methods, a condition on the updated sensitivity along with curve fitting through polynomial regression. Ho wev er, it should be noted that the con ventional KF algorithm’ s con vergence can only be guaranteed in an ideal linear -Gaussian en vironment. Reference [13] has addressed the marginalized PF design in order to estimate the model parameters in the case of hypertension under SNP administration. The method allows considering linear and nonlinear states to be estimated sep- arately to reduce the computational burden. The present work utilizes a multiple-model square root CKF (MMSRCKF) method to effecti vely estimate the non- linear MAP model parameters in the case of hypotension un- der a v asopressor drug injection. The patient’ s MAP response model is described by a first-order parameter-v arying model with a varying input delay . The multiple-model part ad- dresses the hypothesis testing and the estimation of the input delay . The square root (SR) algorithm employs the Cholesky factorization of the error cov ariance matrix to guarantee its positiv e definiteness during numerical operations [14]. For the verification of the proposed method, data from animal experiments is collected at the Univ ersity of T exas, Medical Branch (UTMB) at Galveston. The estimation results are compared to that of MMEKF reported in [15]. The mathematical notation to be used in the paper is as follows. t denotes the continuous-time domain, and k stands for the discrete-time v ariable. For a stochastic process, x k , E [ x k ] denotes its expected value and N { x k ; b x k | k , P k | k } represents a normal Gaussian probability distribution with the mean of b x k | k and the covariance of P k | k . This paper is organized as follo ws. The MMSRCKF online estimation algorithm is developed in Section II. Section III introduces a first-order time-delayed dynamic model to characterize the MAP response to the drug infusion. Section IV presents the estimation results and ev aluation of the per- formance of the proposed MMSRCKF method in comparison to MMEKF . Final remarks are provided in Section V. I I . E S T I M A T I O N P R E L I M I N A R I E S A N D M E T H O D O L O G Y In this section, first, a deri vati ve-free on-line sequential state estimator kno wn as the square root CKF (SRCKF) algorithm is formulated for a general nonlinear discrete-time stochastic system. Subsequently , we formulate the multiple- model (MM) approach and couple it with the introduced SRCKF algorithm for the time delay estimation of a system with an input delay . A. SRCKF Algorithm The Bayesian-based CKF scheme aims at estimating the states of a dynamical system using a probabilistic framew ork [6]. The original CKF state estimation process is susceptible to numerical problems such as indefinite error covariance ma- trix, di vergence phenomenon, and filter instability . In order to tackle these obstacles, we enhance CKF with the square root computation, i.e. the cov ariance matrix is decomposed using a factorization method, such as the Cholesky factorization to guarantee positiv e definiteness within numerical operations [16]. The resulting square roots of the error covariance matrices propagate through the sequential state estimation process. Next, the third-degree spherical-radial rule is used to approximate the multidimensional inte grals inv olved in the Bayesian filtering [17]. Let us consider a general nonlinear discrete-time stochastic system as follows  x k +1 = f ( x k , u k ) + w k , y k = h ( x k , u k ) + v k , k = 0 , 1 , . . . , k f , (1) where x k ∈ R n stands for the unmeasured state vector of the system, u k ∈ R n u is the input vector , and y k ∈ R n y is the measurement vector at the time k , and k f is the final time. f ( x k , u k ) : ( R n , R n u ) 7→ R n and h ( x k , u k ) : ( R n , R n u ) 7→ R n y are known general nonlinear vector mappings, and w k ∈ R n and v k ∈ R n y are statistically independent zero-mean Gaussian process and measurement noise signals, respectiv ely . The probability distribution functions (PDFs) of the noise vectors, namely p ( w k ) and p ( v k ) are kno wn, as well as, the initial state vector PDF p ( x 0 ) . SRCKF seeks to find the estimation of the state vec- tor in the form of a conditional PDF , p ( x k | y k ) , that has the entire knowledge about the current state vector , x k , giv en the entire measurement vectors sequence, i.e. y k = [ y 0 y 1 . . . y k ] . Howe ver , in some cases, a Gaussian approximation of the conditional PDF allows to only com- pute the first two conditional moments, i.e. the mean b x k | k = E [ x k | y k ] and the error co variance matrix P k | k = cov [ x k | y k ] which results in p ( x k | y k ) ≈ N { x k ; b x k | k , P k | k } . By assuming Gaussian white noise vectors, the prediction step (state prediction) and correction step (measurement update) are carried out via integrating a nonlinear function concerning a normal distribution, or b x k +1 | k = E [ x k +1 | y k ] = Z R n f ( x k , u k ) p ( x k | y k ) d x k ≈ Z R n f ( x k , u k ) N { x k ; b x k | k , P k | k } d x k , (2) b y k +1 | k = E [ y k +1 | x k +1 ] = Z R n h ( x k +1 , u k +1 ) p ( y k +1 | x k +1 ) d x k +1 ≈ Z R n h ( x k +1 , u k +1 ) N { x k +1 ; b x k +1 | k , P k +1 | k } d x k +1 . (3) The third-degree spherical-radial rule is utilized to compute the numerical approximation of the moment integrals (2) and (3). Next, for an arbitrary function g ( x ) with Σ as the cov ariance of x , the integral I ( g ) = √ 2 π | Σ | − 1 2 Z R n g ( x ) exp  − 1 2 ( x − µ ) T Σ − 1 ( x − µ )  d x , (4) in the spherical coordinate system becomes I ( g ) = (2 π ) − n 2 Z ∞ r =0 Z U n g ( C r z + µ ) d z r n − 1 e − r 2 2 d r , (5) where x = C r z + µ with k z k = 1 , µ is the mean and C is the Cholesky factor of the cov ariance, Σ , and U n is the unit sphere. Then, we used the symmetric spherical cubature rule to further approximate the integral as I ( g ) = 1 2 n 2 n X i =0 g ( √ n ( C ξ i + µ )) , (6) where ξ i denotes the i th cubature point at the intersection of the unit sphere and its axes. The main benefit of this scheme is that the cubature points are obtained off-line using a third- degree cubature rule [18]. W e follow the steps introduced next to compute the estimation of the state vector via the SRCKF algorithm. SRCKF algorithm: 1) Initialization : The state initial condition is giv en by x 0 | 0 ≡ x 0 with b x 0 = E [ x 0 ] where the initial cov ariance matrix is P 0 | 0 . W e decompose it as P 0 | 0 = S 0 | 0 S T 0 | 0 through the Cholesky factorization, i.e. S 0 | 0 = chol { [ x 0 − b x 0 ][ x 0 − b x 0 ] T } . Then, generate the cubature points, ξ i , for the initia state vector and the fixed weights, w i = w = 1 2 n , for i = 1 , 2 , . . . , 2 n . 2) Time update (Prediction) ( k = 1 , 2 , . . . , k f ) : a) Evaluation of the cubature points X i,k − 1 | k − 1 = S k − 1 | k − 1 ξ i + b x k − 1 | k − 1 . (7) b) Evaluation of the propagated cubature points via the system dynamics X ∗ i,k | k − 1 = f k ( X i,k − 1 | k − 1 , u k − 1 ) . (8) c) Evaluation of the predicted states based on the gen- erated weights and propagated points b x k | k − 1 = 2 n X i =1 w i X ∗ i,k | k − 1 . (9) d) Evaluation of the square root of the covariance matrix of the predicted state error covariance S k | k − 1 = tr iangl e  [ χ ∗ k | k − 1 , S Q k − 1 ]  , (10) where χ ∗ k | k − 1 is a centered weighted matrix, i.e. χ ∗ k | k − 1 = 1 √ 2 n [ X ∗ 1 ,k | k − 1 − b x k | k − 1 X ∗ 2 ,k | k − 1 − b x k | k − 1 · · · X ∗ 2 n,k | k − 1 − b x k | k − 1 ] , (11) and S Q k − 1 is the square-root of the the process noise such that Q k − 1 = S Q k − 1 S T Q k − 1 . Moreover , B = tr iangl e { A } stands for a general triangularization algorithm, e .g. QR decomposition, where B is a lower triangular matrix. If C is an upper triangular matrix obtained through the QR decomposition of A T , then the lower triangular matrix is gi ven by B = C T . 3) Measurement update (Correction) ( k = 1 , 2 , . . . , k f ) : a) Evaluation of the cubature points using the predicted square root matrix, S k | k − 1 , X i,k | k − 1 = S k | k − 1 ξ i + b x k | k − 1 . (12) b) Evaluation of the propagated cubature point via the output dynamics Y i,k | k − 1 = h ( X i,k | k − 1 , u k ) . (13) c) Estimation of the predicted measurement vector b y k | k − 1 = 2 n X i =1 w i Y i,k | k − 1 . (14) d) Evaluation of the square root of the innovation co- variance matrix S y y,k | k − 1 = tr iangl e  [ Y k | k − 1 , S R k ]  , (15) where Y k | k − 1 is a centered weighted matrix, i.e. Y k | k − 1 = 1 √ 2 n [ Y 1 ,k | k − 1 − b y k | k − 1 Y 2 ,k | k − 1 − b y k | k − 1 · · · Y 2 n,k | k − 1 − b y k | k − 1 ] . (16) S R k is also the square-root of the the measurement noise such that R k = S R k S T R k . e) Evaluation of the cross-cov ariance matrix P xy ,k | k − 1 = χ k | k − 1 Y T k | k − 1 , (17) with the centered weighted matrix χ k | k − 1 giv en by χ k | k − 1 = 1 √ 2 n [ X 1 ,k | k − 1 − b x k | k − 1 X 2 ,k | k − 1 − b x k | k − 1 · · · X 2 n,k | k − 1 − b x k | k − 1 ] . (18) f) Evaluation of the SRCKF filter gain W k = P xy ,k | k − 1 S − T y y,k | k − 1 S − 1 y y,k | k − 1 . (19) g) Evaluation of the corrected state update based on the measurement b x k | k = b x k | k − 1 + W k ( y k − b y k | k − 1 ) . (20) h) Evaluation of the square-root of the corrected error cov ariance matrix S k | k = tr iangl e  [ χ k | k − 1 − W k Y k | k − 1 , W k S R k ]  . (21) The state estimation process continues iterativ ely from the second step of the algorithm, i.e. the time update (prediction) by setting k = k + 1 . B. Multiple-Model SRCKF for Input Delay Estimation T ime delay estimation introduces a challenge in the param- eter identification frame work since the variable delay is not transformable to an equiv alent random walk process. Ratio- nal approximations of the delay such as Pad ´ e approximation can be considered as alternativ e solutions; howe ver , the introduced truncation error may be significant and problem- atic, especially for large and time-v arying delays. Thus, to obtain a more accurate delay estimation, the aforementioned SRCKF algorithm is equipped with a multiple-model (MM) framew ork cascaded with a hypothesis testing module [19]. The underlying idea of the MMSRCKF method is to use a bank of N identical SRCKFs in a parallel setting, as shown in Fig. 1. Every SRCKF uses the same measurement and input data, but a dif ferent delay is assigned to each filter . . . . . . . SRCKF with τ i SRCKF with τ 1 SRCKF with τ N Hypothesis T esting u k y k b X 1 k b X i k b X N k b τ M M k Fig. 1: Bank of N parallel SRCKFs for delay estimation The i th element in the bank provides us with a state vector estimation X i k together with the residuals r i k = y k − b y i k . By having this information, a hypothesis testing block can then be used to estimate the value of the delay . Specifically , if the delay matches the one assigned to the i th SRCKF element, then the corresponding residual is essentially a zero-mean white noise process, i.e. E [ r i k ] = 0 , and its cov ariance is giv en by E [ r i k ( r i k ) T ] = HP i k H T + R , R i k , (22) where H = [1 0 0 1] , P i k denotes the estimation cov ariance at the k th step, and R denotes the measurement noise co- variance. The conditional probability density function of the i th SRCKF element measurement can be computed through f ( b y i k | y k ) = 1 (2 π ) m 2 | R i k | 1 2 exp n − 1 2 ( r i k ) T ( R i k ) − 1 r i k o , (23) where m is the dimension of av ailable measurements at each time step. Then, the conditional probability of each hypothesis is p i k = f ( b y i k | y k ) p i k − 1 N P j =1 f ( b y j k | y k ) p j k − 1 , (24) where p i k can be interpreted as the normalized conditional probability of the case when the delay equals the assigned value to the i th filter, i.e. N P j =1 p j k = 1 . Now , it is possible to estimate the delay according to the filter , which has the highest probability . Howe ver , to obtain a more accurate delay estimation and to a void large fluctuations, instead of choosing the block with most likely delay estimation, we treat the hypotheses resulting as weights and blend them to improve the delay estimation. In other words, we can estimate the time delay as ˆ τ M M k = N X j =1 p j k τ j k , (25) where τ j k is the delay estimation of the i th filter . In the next section, we will present the mathematical model describing the dynamics of the MAP response to the PHP drug infusion. Fig. 2: T ypical MAP variations in response to PHP step injection [15] I I I . M A P R E S P O N S E M O D E L I N G The following first-order model with an input delay has been broadly considered and implemented in the literature to characterize the patient’ s MAP response to the infusion of a vasoacti ve drug, such as phenylephrine (PHP) [4], [15], [20]: T ( t ) · ˙ ∆ M AP ( t ) + ∆ M AP ( t ) = K ( t ) · u ( t − τ ( t )) , (26) where ∆ M AP ( t ) stands for the MAP changes in mmH g from its baseline value, i.e. ∆ M AP ( t ) = M AP ( t ) − M AP b ( t ) , u ( t ) is the drug injection rate in ml/h , K ( t ) denotes the patient’ s sensitivity to the administered drug, T ( t ) is the lag time describing the uptake, distribution and biotransformation of the drug [21], and τ ( t ) is the time delay for the drug to reach the circulatory system from the injection site. This model structure seems to adequately describe a patient’ s physiological response to the PHP drug injection. Fig. 2 presents a typical MAP response due to a step PHP infusion v ersus a matched response of (26). This figure also illustrates the interpretation of the model parameters K ( t ) , T ( t ) , τ ( t ) , M AP b ( t ) which hav e been obtained to fit the MAP response using a least-squares optimization method. Although the proposed model structure (26) is qualitativ ely able to represent the characteristics of the MAP response to the infusion of PHP , experiments show that the model parameters vary significantly over time due to patients’ pharmacological variability subject to the vasoacti ve drug infusion. That is, the model parameters and delay could vary remarkably from patient-to-patient (inter- patient variability), as well as, for a given patient ov er time (intra-patient variability) [21], [22]. For the implementation of recursiv e sequential estimation tools, we discretize the continuous-time model (26) at the sampling rate of T s as follows    x k +1 =  1 − T s T k  x k + K k T s T k u ( k − τ k T s ) , y k = x k + M AP b k , (27) where x k = ∆ M AP k = M AP k − M AP b k at the k th time instant. In (27), we augment the state vector with the parameters to be estimated, namely K k , T k , and M AP b k by assuming local random-walk dynamics. In other words X T k = [ X 1 k X 2 k X 3 k X 4 k ] = [∆ M AP k K k T k M AP b k ] . (28) Since model parameters are all time-varying and assumed to be a priori unknown, (27) represents a nonlinear equation with regards to the state vector , X k , that can be expressed as the following nonlinear dynamics ( X 1 k +1 = f k ( X k , u k ) + w k , y k = h k ( X k ) + v k , (29) with    f 1 k ( X k , u k ) =  1 − T s X 3 k  X 1 k + T s X 2 k X 3 k u ( k − τ k T s ) , h k ( X k ) = X 1 k + X 4 k . (30) The process noise, w k , and the measurement noise, v k , are both assumed to be additiv e and statistically independent zero-mean Gaussian processes with cov ariances given by Q k and R k , respecti vely . Although such an augmentation facili- tates the estimation procedure, the time-varying input delay neither can be included in the augmented state vector nor be captured by a random walk process. Thus, the time delay is estimated through a multiple-model (MM) hypothesis testing process along with the SRCKF , discussed in Section II. Next, we test the proposed MMSRCKF estimation algo- rithm in a simulation where the patient’ s model parameters are generated by nonlinear functions based on clinical obser- vations. Then, we validate the verified estimation framework using the experimental data from animal experiments. I V . M A P R E S P O N S E E S T I M A T I O N R E S U L T S A N D V A L I D A T I O N S In order to v alidate the proposed parameter estimation method, first, we need to b uild a realistic simulation model of an individual’ s MAP response to the drug infusion with some known model parameters. Adopting (26), we generate the patient’ s nonlinear time-varying model parameters, i.e. K ( t ) , T ( t ) , τ ( t ) , and M AP b ( t ) based on clinical observations as follows [4] a k ˙ K ( t ) + K ( t ) = k 0 exp {− k 1 i ( t ) } , T ( t ) = sat [ T min ,T max ] { b T Z t 0 i ( t ) dt } , ( a τ , 2 ... τ ( t ) + a τ , 1 ¨ τ ( t ) + ˙ τ ( t ) = b τ , 1 ˙ i ( t ) + i ( t ) , t ≥ t i 0 , τ ( t ) = 0 , other wise, (31) where i ( t ) is the drug injection and a k , k 0 , k 1 , b T , a τ , 2 , a τ , 1 , and b τ , 1 are uniformly distributed random coef ficients gi ven in T able I [1]. Also, the MAP baseline value, M AP b ( t ) , is assumed to be constant and equal to 70 mmH g for the con- sidered nonlinear patient. As per (31), the model parameters K ( t ) , T ( t ) , and τ ( t ) are nonlinear functions of the drug infusion rate, i ( t ) . Fig. 3 demonstrates the general structure of the nonlinear patient parameter generation process. The model parameters are generated based on the giv en infusion T ABLE I: Probabilistic distributions of coef ficients in (31) Parameter Distribution a k U (500 , 600) k 0 U (0 . 1 , 1) k 1 U (0 . 002 , 0 . 007) b T U (10 − 4 , 3 × 10 − 4 ) a τ , 1 U (5 , 15) a τ , 2 U (5 , 15) b τ , 1 U (80 , 120) Fig. 3: Structure of nonlinear patient parameter generation T ABLE II: Estimation root mean square errors (RMSEs) RMSE Parameter MMSRCKF EKF K 0.061 0.095 T 8.370 24.917 M AP b 0.188 0.706 τ 3.128 9.114 M AP 0.008 0.202 rate, i ( t ) , while the parameter estimation tool estimates the model parameters sub-optimally , using the input drug infusion rate and measured output MAP . Fig. 4 shows the piecewise constant PHP drug infusion profile, that is used to generate the nonlinear patient parameters. Using the generated model parameters, we e valuate the performance of the proposed MMSRCKF method in estimating them, and the estimation results are compared to the previously reported EKF algorithm [15]. Figs. 5, 6, 7, and 8 show the estimation results for the model parameters, namely the sensitivity K ( t ) , time constant T ( t ) , MAP baseline value M AP b ( t ) , and time delay τ ( t ) , respectiv ely . As we can see, the implemented MMSRCKF method outperforms the EKF in terms of accuracy and the con vergence speed. The MMSR- CKF online estimation results show better matches with the generated nonlinear patient reference parameters. It should be noted that the computation complexity of both CKF and EKF algorithms equally grows as n 3 where n denotes the system size where the former filter is more accurate and numerically more stable. T able II further compares the root mean square errors (RMSEs) of the model parameters and estimated MAP response in both algorithms by which the error reduction is obvious using MMSRCKF . Fig. 4: Profile of piecewise constant PHP drug injection Fig. 5: Nonlinear patient sensitivity estimation Fig. 6: Nonlinear patient lag time estimation In the next step, we implemented the MMSRCKF al- gorithm on the collected data from an actual animal ex- periment. The input PHP drug infusion rates and output MAP measurements for a 55 kg anesthetized swine were recorded at the Resuscitation Research Laboratory at the Department of Anesthesiology , UTMB at Galveston, T exas. Precisely speaking, an intramuscular injection of ketamine was used to sedate the swine which were maintained under anesthetic conditions by the continuous infusion of propofol. A Philips MP2 transport device with a sampling frequency of 20 H z was used to monitor the blood pressure response ov er a 6 -hour experiment, while the PHP drug was being infused through a bodyguard infusion pump. Fig. 9 shows Fig. 7: Nonlinear patient baseline MAP estimation Fig. 8: Nonlinear patient input delay estimation the piecewise constant PHP drug infusion profile versus the corresponding measured raw blood pressure response and the MAP response over time. W e then utilized this dataset for the validation of the estimation of the MAP dynamic model parameters using the proposed MMSRCKF methodology . The experimental dataset was re-sampled at the sampling frequency of 0 . 2 H z . Regarding the multiple-model part of the MMSRCKF algorithm for the delay estimation, the esti- mation accuracy versus the algorithm speed of con ver gence triggered a trade-off which needed to be addressed with care; hence, it was essential to choose an appropriate number of the bank of SRCKFs constructing the MMSRCKF structure. In this work, we examined a bank of 11 SRCKFs with the delay interv al of τ ∈ [0 100] s . Consequently , the time gridding for the ev enly distributed filters was equal to 10 s . The MAP estimation of the proposed MMSRCKF algorithm, as well as the clinically acquired MAP measurements, are illustrated in Fig. 10, which suggests that the proposed identification method is capable of accurately capturing the MAP response of the swine to the injection of the PHP drug. Additionally , the estimation of the model parameters, namely the sensitivity K ( t ) , time constant T ( t ) , MAP baseline value M AP b ( t ) , and time delay τ ( t ) , are depicted in Figs. 11, 12, 13, and 14, respectiv ely . The estimated parameter values followed the expected trends, as discussed in detail in [4]. Furthermore, the delay estimation in Fig. 14 demonstrated a sharp initialization peak right after the initial injection of the Fig. 9: Instantaneous blood pressure and MAP response to a piecewise constant PHP drug injection in animal experiments Fig. 10: MAP estimation results in animal experiments Fig. 11: Sensitivity estimation in animal experiments drug and follo wed a slowly decaying trend during the rest of the experiment as anticipated [1]. V . C O N C L U S I O N Precise estimation of hemodynamics characteristics and mean arterial blood pressure in response to vasoacti ve drug administration is pi votal to design an effecti ve controller to meet closed-loop physiological response requirements in various clinical scenarios. Real-time estimation of such dynamic models was e xamined in this paper . Due to the inter- and intra-patient variability , a parameter-v arying model with Fig. 12: Lag time estimation in animal experiments Fig. 13: Baseline MAP estimation in animal experiments Fig. 14: Time delay estimation in animal experiments varying input delay was deemed to account for model pa- rameter v ariations. A Bayesian estimation scheme known as cubature Kalman filter was used because of its con vergence speed, nonlinear system handling, and numerical stability . The varying parameters of the nonlinear system corrupted by noise were estimated through the proposed framework. Since the input delay cannot be captured via a random-walk pro- cess, the filter was augmented with a multiple-model module. Delay and parameters estimation results in comparison to classical extended Kalman filter were demonstrated which verified the advantage of the utilized Bayesian approach. R E F E R E N C E S [1] C. R. Craig and R. E. Stitzel, Modern pharmacology with clinical applications . Baltimore, MD, USA: Lippincott Williams & W ilkins, 2004. [2] S. J. da Silva, T . A. Scardovelli, S. R. M. da Silva Boschi, S. C. M. Rodrigues, and A. P . da Silva, “Simple adapti ve PI controller develop- ment and evaluation for mean arterial pressure regulation, ” Research on Biomedical Engineering , vol. 35, no. 2, pp. 157–165, 2019. [3] S. Herget-Rosenthal, F . Saner , and L. S. Chawla, “ Approach to hemo- dynamic shock and vasopressors, ” Clinical Journal of the American Society of Nephr ology , vol. 3, no. 2, pp. 546–553, 2008. [4] S. T asoujian, S. Salavati, M. Franchek, and K. Grigoriadis, “Robust imc-pid and parameter-v arying control strategies for automated blood pressure regulation, ” International Journal of Control, Automation and Systems , vol. 17, no. 7, pp. 1803–1813, 2019. [5] K. Kashihara, T . Kaw ada, K. Uemura, M. Sugimachi, and K. Suna- gawa, “ Adaptive predictiv e control of arterial blood pressure based on a neural network during acute hypotension, ” Annals of Biomedical Engineering , vol. 32, no. 10, pp. 1365–1383, 2004. [6] I. Arasaratnam and S. Haykin, “Cubature Kalman filters, ” IEEE T ransactions on Automatic Contr ol , vol. 54, no. 6, pp. 1254–1269, 2009. [7] X. Zhao, J. Li, X. Y an, and S. Ji, “Robust adaptive cubature Kalman filter and its application to ultra-tightly coupled SINS/GPS navigation system, ” Sensors , vol. 18, no. 7:2352, pp. 1–19, 2018. [8] E. Furutani, M. Araki, S. Kan, T . Aung, H. Onodera, M. Imamura, G. Shirakami, and S. Maetani, “ An automatic control system of the blood pressure of patients under surgical operation, ” International Journal of Contr ol, Automation, and Systems , vol. 2, no. 1, pp. 39–54, 2004. [9] Y . Gao and M. J. Er , “ An intelligent adaptiv e control scheme for postsurgical blood pressure regulation, ” IEEE T ransactions on Neural Networks , vol. 16, no. 2, pp. 475–483, 2005. [10] J. B. Slate, L. C. Sheppard, V . C. Rideout, and E. H. Blackstone, “ A model for design of a blood pressure controller for hypertensiv e patients, ” IF AC Pr oceedings , vol. 12, no. 8, pp. 867–874, 1979. [11] K. Y . Zhu, H. Zheng, and D. G. Zhang, “ A computerized drug delivery control system for regulation of blood pressure, ” International Journal of Intelligent Computing in Medical Sciences & Image Processing , vol. 2, no. 1, pp. 1–13, 2008. [12] N. Malagutti, A. Dehghani, and R. A. Kennedy , “Robust control design for automatic regulation of blood pressure, ” IET Control Theory & Applications , vol. 7, no. 3, pp. 387–396, 2013. [13] N. Malagutti, “Particle filter-based robust adaptiv e control for closed- loop administration of sodium nitroprusside, ” Journal of Computa- tional Surgery , v ol. 1, no. 1, pp. 1–19, 2014. [14] X. Cui, Z. He, E. Li, A. Cheng, M. Luo, and Y . Guo, “State-of-char ge estimation of power lithium-ion batteries based on an embedded micro control unit using a square root cubature Kalman filter at various ambient temperatures, ” International J ournal of Ener gy Resear ch , vol. 43, no. 8, pp. 3561–3577, 2019. [15] T . Luspay and K. M. Grigoriadis, “ Adaptiv e parameter estimation of blood pressure dynamics subject to v asoactiv e drug infusion, ” IEEE T ransactions on Contr ol Systems T echnology , vol. 24, no. 3, pp. 779– 787, 2015. [16] N. Loehr , Advanced linear algebra . New Y ork, NY , USA: Chapman and Hall/CRC, 2014. [17] B. Jia, M. Xin, and Y . Cheng, “High-degree cubature Kalman filter, ” Automatica , vol. 49, no. 2, pp. 510–518, 2013. [18] Y . Liu, K. Dong, H. W ang, J. Liu, Y . He, and L. Pan, “Adaptive Gaussian sum squared-root cubature Kalman filter with split-merge scheme for state estimation, ” Chinese Journal of Aer onautics , v ol. 27, no. 5, pp. 1242–1250, 2014. [19] P . D. Hanlon and P . S. Maybeck, “Multiple-model adaptive estimation using a residual correlation Kalman filter bank, ” IEEE T ransactions on Aer ospace and Electronic Systems , vol. 36, no. 2, pp. 393–406, 2000. [20] C. Sandu and D. Popescu, “Reinforcement learning for the control of blood pressure in post cardiac surgery patients, ” U. P . B. Sci. Bull., Series C , vol. 78, no. 1, pp. 139–150, 2016. [21] S. Isaka and A. V . Sebald, “Control strategies for arterial blood pressure regulation, ” IEEE T ransactions on Biomedical Engineering , vol. 40, no. 4, pp. 353–363, 1993. [22] R. R. Rao, B. Aufderheide, and B. W . Bequette, “Experimental studies on multiple-model predicti ve control for automated re gulation of hemodynamic variables, ” IEEE T ransactions on Biomedical Engi- neering , vol. 50, no. 3, pp. 277–288, 2003.