From the *Research Division on Advanced Information Technology, Information Synergy Center, Tohoku University, Sendai, Japan; †Institute of Engineering Mechanics and Systems, University of Tsukuba, Tsukuba, Japan; ‡Department of Electrical and Communication Engineering, Graduate School of Engineering, Tohoku University, Sendai, Japan; §Department of Electrical Engineering, Tohoku Gakuin University, Tagajo, Japan; ¶Institute of Development, Aging and Cancer, Tohoku University, Sendai, Japan; and ∥Department of Surgery, Baylor College of Medicine, Houston, Texas.

Submitted for consideration July 2000; accepted for publication in revised form January 2002.

Reprint requests: Makoto Yoshizawa, Research Division on Advanced Information Technology, Information Synergy Center, Tohoku University, Electrical and Information Building, Aoba-yama 05, Aoba-ku, Sendai 980-8579, Japan.

For automatic control of an implantable artificial heart, it is desirable to estimate blood flow and pressure without using flowmeters and pressure sensors implanted inside the recipient’s body. In the case of an electric motor driven artificial heart, voltage, current, and rotational speed of the motor can be measured easily and accurately. Because these three signals are closely related to flow of the blood pump and the pressure difference (pressure head) between inlet and outlet ports, it is possible to estimate flow and pressure head by processing these signals. Such an approach has already succeeded to some extent in the case of centrifugal blood pumps in other studies. ^{1–5} The methods of estimation based on this approach must include a compensation procedure for the change in blood viscosity. However, conventional compensation procedures are not suitable for automatic and on-line estimation. For example, hematocrit value may be useful to compensate for the change in blood viscosity, ^{3–5} but direct measurement of this value should be avoided because of the necessity for another measurement apparatus outside the body.

The present study has proposed a new method for estimating instantaneous values of flow (Q (t )[L/min]) and pressure head (P (t )[mm Hg]) on the basis of voltage (V (t )[V]), current (I (t )[A]), and rotational speed (N (t )[k(rpm)]) in a DC motor driven centrifugal pump used as a ventricular assist device (VAD). The proposed estimator can automatically compensate for the effect of the change in blood viscosity on the estimation accuracy by employing two auto-regressive exogenous (ARX) models. ^{6}

Methods
Steady State Model
Define electrical power (VI (t )[W]) supplied to the motor as VI (t ) =V (t )·I (t ). As shown in Figure 1 , the relationship between flow Q and power VI parameterized by rotational speed N could be plotted on the basis of the data obtained from an in vitro experiment that will be described in detail later. This figure means that, in a static situation, the relationship can be approximated by a line if N is kept constant, and that the slope (α_{flow} (N )) and the bias (the flow-axis intercept; β_{flow} (N )) of the line can be represented as functions of N . In fact, the slope and the bias could be approximated by the second order functions of N as shown in Figure 2 . It has already been ascertained that the similar approximation can also hold in the relationship among pressure head P , power VI , and rotational speed N using the slope (α_{pressure} (N )) and the bias (β_{pressure} (N )).

Figure 1: Relationship between flowQ and power VI parameterized by rotational speed N plotted on the basis of the data obtained from an in vitro experiment.

Figure 2: The slope and bias of the approximate line shown inFigure 1 plotted versus the rotational speed N. The plots can be approximated by the second order polynomial functions of N.

For general expression, let y denote either P or Q , and then the static approximation described above can be expressed bywhere the slope α(N ) and the bias β(N ) are second order polynomial functions of N whose subscript “flow” or “pressure” is omitted here:where b _{j} ;j = 1, 2, ...,6 is a constant to be identified by measured data. Expanding equation 1 using equations 2 and 3 , we have

ARX Model I for Output Estimation without Additional Input
To estimate the transient response of flow or pressure head, the static system (equation 1 ) should be extended to the following ARX model (ARX model I ) as a dynamic time series model illustrated in Figure 3A . where k is the discrete time satisfying t =kΔt , Δt is the sampling period, w (k ) is the residue assumed to be a white noise, n is the order of the output, m _{j} is the order of the input, and u _{j} (k ) is six kinds of exogenous input corresponding to equation 4 as follows:Let m _{6} = 1 and m _{1} =m _{2} = ... =m _{5} =m for simplicity.

Figure 3: Two auto-regressive exogenous (ARX) models to estimate flowQ (t) or pressure head P (t). A: ARX model I estimates on the basis of supplied power VI (t) and rotational speed N (t). B: The combination of ARX models II and III estimates on the basis of VI (t), N (t), and the gain constant K (t), in which ARX model III yields K (t) on the basis of VI (t) and N (t). K (t) is related to rotational speed per unit power and reflects the change in viscosity of blood or hematocrit.

Before implantation of an artificial heart, we can collect the measured data y (k ) and u _{j} (k ). Thus, the coefficient parameters a _{i} and b _{ij} included in equation 5 can be identified by an appropriate batch (off-line) identification method ^{6} on the basis of y (k ) and u _{j} (k ). The off-line least squares method ^{6} was used in the present study. After the identification of a _{i} and b _{ij} , the ARX model I given by equation 5 without w (k ) can generate an estimate (ŷ (k )) of y (k ) on the basis of on-line measurement of u _{j} (k ).

In general, the final prediction error ^{6} or the Akaike’s information criterion ^{6} is frequently used to determine the orders (n and m ) of ARX models. In this case, however, it has already been ascertained that the determination based on such a criterion tends to yield extremely high orders. A system model with too high of orders is apt to produce a large estimation error, even if the true system slightly changes from the original identified situation. To avoid such an over-learning phenomenon, the number n and m can be chosen as under 20 so as to minimize the root mean square value (ε) of the error between y (k ) and ŷ (k ) defined aswhere K _{D} is the size of data.

In the present study, ε was used to evaluate the estimation accuracy of the model. Moreover, the correlation coefficient r defined by equation 13 was used to evaluate linear correlation between y (k ) and ŷ (k ).where JOURNAL/asaio/04.02/00002480-200207000-00020/ENTITY_OV0374/v/2017-07-29T043104Z/r/image-pngy and ŷ are mean values of y (k ) and ŷ (k ), respectively.

ARX Model II for Output Estimation with Additional Input
If blood viscosity changes considerably, the estimation accuracy of the ARX model I given by equation 5 will decrease because the true coefficients a _{i} and b _{ij} will also change. Thus, as shown in the right block in Figure 3B , let us consider a new ARX model (ARX model II ) as follows:where the first six inputs are the same as equations 6 to 11 , respectively, and the seventh input is a new additional input given byMATH where K is the steady state gain obtained from ARX model III (the left block in Figure 3B ), which will be explained later. Let m _{6} =m _{7} = 1 and m _{1} =m _{2} = ... =m _{5} =m . The method of determining m is the same as that of ARX model I.

ARX Model III for Estimation of Steady State Gain
The change in blood viscosity must affect a parameter associated with the friction force of the motor-pump system. Therefore, we can use such a parameter as an additional input to compensate for the change in the true system model because the parameter includes information on blood viscosity that may prevent rotational motion. One candidate for such a parameter is the steady state gain K of the system from the supplied power VI (t ) to the rotational speed N (t ). The gain K is directly related to rotational speed per unit torque associated with the given electrical power, and K then correlates with the reciprocal of the friction coefficient. Assume that the system from VI (t ) to N (t ) can be approximated by another very simple ARX model (ARX model III ):Let b _{12} = 0, w (k ) = 0. For k → ∞, we can regard N (∞) =N (k ) =N (k - 1) =N (k - 2) and VI (∞) =VI (k ) =VI (k - 1); the steady state gain K =N (∞)/ VI (∞) can be identified by using the coefficients of equation 16 as follows:MATH

This identification process can be carried out every time the viscosity seems to change because both VI (k ) and N (k ) are always measurable inside the body.

In this study, the function arx in the System Identification Toolbox of Matlab (Mathworks, Inc., Natick, MA ) was used to identify the coefficients of ARX models I, II, and III.

Experiment
As shown in Figure 4 , an in vitro experiment was carried out in a mock circulatory system consisting of a centrifugal pump (Kyocera C1E3; Kyoto, Japan) driven by a brushless DC motor, tubes, a reservoir, and a clamp. To compare different effects of change in viscosity on the estimation accuracy, two different solutions were used in the mock system, i.e., glycerin 37% and water. An electromagnetic flowmeter was used to measure the outflow of the pump. Two pressure sensors equipped at the inlet and outlet ports of the pump were used to measure pressure head. Sampling period Δt was 10 ms. To maintain the persistent excitation condition ^{6} for accurate identification, the driving voltage V (t ) was randomly changed between 5V and 9V. The gap of the clamp was also randomly changed by hand to simulate the effect of transient change in cardiac load made by the natural heart connected in parallel with the blood pump as the VAD.

Figure 4: Schematic illustration of the mock circulatory system driven by a centrifugal pump.

Results
Figure 5 shows time trajectories of V (t ), I (t ), VI (t ), N (t ), P (t ), and Q (t ). The left and the right halves of Figure 5 correspond to glycerin 37% (interval A) with high viscosity and water (interval B) with low viscosity, respectively. The data file of Figure 5 was made by simply connecting two data files corresponding to intervals A and B that had been independently measured at different times.

Figure 5: Time trajectories of measurementsV (t), I (t), VI (t), N (t), Q (t), and P (t) when the motor was driven by randomly changed voltage V (t) between two levels, and the outlet tube was randomly clamped. Left and right halves correspond to 37% glycerin and water, respectively.

Figure 6 shows time trajectories of the input u _{j} (t ), j = 1, 2, ...,7 except j = 6. The value of the steady state gain K was identified as K = 0.278k(rpm)/W in interval A and K = 0.300k(rpm)/W in interval B. This means that K increased by about 8% given the reduction of viscosity from glycerin to water.

Figure 6: Time trajectories of inputs used by auto-regressive exogenous (ARX) models,N2(t)·VI (t), N (t)·VI (t), VI (t), N2(t), N (t) and K (t) in the same situation as Figure 5.

The dotted line in Figure 7 shows the result of estimation of outflow Q (t ) using ARX model I defined by equation 5 whose coefficients were identified on the basis of the data obtained from only interval A. Thus, it is a matter of course that Q (t ) could be well estimated in the same interval A but considerably underestimated in interval B in which the corresponding input–output model must have different coefficients.

Figure 7: Flow estimation using auto-regressive exogenous (ARX) model I (n = 12,m = 12) whose coefficients were identified in only Interval A. Large underestimation was found in Interval B. Root mean square error: ε = 1.99L/min. Correlation coefficient:r = 0.836.

On the contrary, Q (t ), depicted as the dotted line in Figure 8 , was estimated using ARX model I, whose coefficients were identified on the basis of the data obtained from only interval B. Large overestimation can be observed in interval A, whose data were not used for identification, while the estimated value agrees well with the measured value in interval B.

Figure 8: Flow estimation using auto-regressive exogenous (ARX) model I (n = 13,m = 13) whose coefficients were identified in only Interval B. Large overestimation was found in Interval A. ε = 2.07L/min. r = 0.812.

On the other hand, Figure 9 shows the estimate of Q (t ) using ARX model I whose coefficients were identified on the basis of the data obtained from both intervals A and B. It can be seen that Q (t ) could not be well estimated in both intervals. The error ε defined by equation 12 and the correlation coefficient r defined by equation 13 were 2.22L/min and 0.696, respectively.

Figure 9: Flow estimation using auto-regressive exogenous (ARX) model I (n = 10,m = 10) with coefficients identified in both Intervals A and B. The estimated flow did not agree well with that measured in both Intervals A and B. ε = 2.22L/min. r = 0.696.

Figure 10 shows the estimate of Q (t ) using a combination of ARX model II defined by equation 14 and ARX model III defined by equation 16 . First, ARX model III was identified in each interval to provide the corresponding two estimated values of K . Second, the coefficients of ARX model II were identified on the basis of K s as well as the data obtained from both intervals A and B. Finally, the estimate of Q (t ) was calculated from the coefficients in both intervals A and B. The error ε and correlation coefficient r were 1.66L/min and 0.850, respectively. This means that the estimation accuracy of Figure 10 was considerably improved by about 25%, in comparison with Figure 9 .

Figure 10: Flow estimation using the combination of auto-regressive exogenous (ARX) models II (n = 12,m = 12) and III (n = 2, m = 2) with coefficients identified during both Intervals A and B. The estimation accuracy was improved in comparison with Figure 9. ε = 1.66 L/min. r = 0.850.

Figure 11 shows the estimate of pressure head P (t ) using the combination of ARX models II and III. The error ε and correlation coefficient r were 12.8 mm Hg and 0.906, respectively. The estimate agreed slightly less well with the measured pressure in high pressure parts than in low pressure parts. However, the correlation coefficient was higher than in the case of flow estimation.

Figure 11: Pressure head estimation using the combination of auto-regressive exogenous (ARX) models II (n = 11,m = 15) and III (n = 2, m = 2) with coefficients identified during both Intervals A and B. ε = 12.8 mm Hg. r = 0.906.

In the four cases made by the combination of two estimates and two models, ε and r obtained from the same data are shown in Table 1 . The table indicates that pressure estimation yielded much less difference between ARX model I and the combination of ARX models II and III than did flow estimation. This means that pressure estimation was little affected by the change in fluid viscosity, at least in the centrifugal pump used in the present study.

Table 1: Estimation Error ε and Correlation Coefficientr

Discussion
If the estimation task is carried out only when the viscosity does not change, the estimation error will be smaller. In clinical settings, however, this situation is not realistic because the viscosity is not always constant. In general, hemodynamic parameters of a recipient in the estimation process are different from those of an animal or a mock circulatory system used for identification of the model. Even if identification is feasible in the same recipient, the viscosity must change according to the change in the recipient’s hemodynamic state. Namely, blood viscosity changes considerably depending not only on hematocrit value, but also on flow rate, because blood is a non-Neutonian fluid. ^{7} The difference between a 37% solution glycerin and water introduced in our experimental setting may not be so extreme in comparison with clinical settings because such difference is likely to be caused by the change in hematocrit and/or flow rate.

Tsukiya et al. ^{1,2} proposed an indirect measurement method for flow and pressure differences in a magnetically suspended centrifugal pump. To overcome the effect of the change in fluid viscosity, their method needs a calibration procedure, during which it is necessary to clamp an outflow tube of the pump and shut down blood flow, which should be avoided because of inducing thrombosis formation. On the other hand, in our method, the ARX model III for obtaining the steady state gain K can be identified without clamping the outlet tube after implantation of the pump; random clamping of the outlet tube is needed only before implantation.

Wakisaka et al. ^{3–5} proposed another method of estimating flow and pressure head that includes a compensation process for the change in viscosity. Because their method uses a hematocrit that needs an off-line measurement procedure, it is difficult to extend their method to a fully automatic estimation system. On the other hand, the steady state gain K reflecting fluid viscosity can be identified easily on the basis of voltage, current, and rotational speed. Because these measurements can be derived directly from the motor, the estimation process of our method has the capability of being carried out in an automatic, on-line fashion.

Our method needs at least two different data sets with variation in viscosity to identify the coefficients of ARX model II given by equation 14 , especially b _{17} with respect to the additional input K in a mock circulatory system before implantation of the VAD. Of course, to obtain a more accurate estimate, fine tuning of this scalar value must be obtained during implantation.

As shown in Table 1 , the correlation coefficient of our method was not high in comparison with other estimation methods. ^{1–5} For example, Wakisaka et al. ^{5} reported that the correlation coefficient between the measured flow and its estimate was r = 0.994 in a static situation. The lower estimation accuracy of our method may be caused by the fact that the test data used in the present experiment were produced under random change in voltage and random clamping at the outlet tube to ascertain estimation ability in the transient state. This task was required because transient behavior will be more important than static behavior if estimated, instead of measured, values are used for controlling the artificial heart as rapidly as possible or for detecting abnormalities such as suction. On the other hand, other estimation methods were evaluated only in the steady state or near steady state. Thus, it should be noted that our method can not be directly compared with other methods in estimation accuracy.

Conclusion
The present study has proposed a new method for estimating flow and pressure head of a centrifugal pump by using voltage, current, and rotational speed of the motor for the pump without any additional sensors. The proposed estimator consists of two ARX models. One ARX model has one output (pressure head or flow) and three inputs (supplied power, i.e., product of voltage and current; rotational speed; and the steady state gain of the system from supplied power to rotational speed). The steady state gain is felt to reflect viscosity of blood and is identified by means of the other ARX model. The proposed method has attained about 25% improvement in estimation accuracy in comparison with another ARX model without the steady state gain as an input. On the other hand, it has been shown that there was little effect of additional input with respect to the steady state gain on the estimation accuracy in the case of pressure estimation.

To improve the estimation accuracy, artificial neural networks ^{8,9} could be introduced to the estimator to compensate for the nonlinear part of the residue that cannot be expressed by linear models, such as ARX, in further studies.

Acknowledgment
The authors thank Mr. Don Flowers, Ms. Julie Glueck, Ms. Mary Jackson, Dr. Shinji Kawahito, Dr. Kenji Nonaka, Dr. Tamaki Takano, and Dr. Tomohiro Maeda for helpful assistance, fruitful discussion, and excellent cooperation. This work has been supported by the Program for Promotion of Fundamental Studies in Health Science of the Organization for Pharmaceutical Safety and Research of Japan and Grant-in-Aid (No. 08558091) for Scientific Research from the Japan Ministry of Education, Science, Culture and Sports.

References
1. Tsukiya T, Akamatsu T, Nishimura K: Flow rate measurement for magnetically suspended centrifugal blood pump. Jpn J Artif Organs 26: 98–102, 1997.

2. Tsukiya T, Akamatsu T, Nishimura K, Yamada T, Nakazeki T: Use of motor current in flow rate measurement for the magnetically suspended centrifugal blood pump. Artif Organs 21: 396–401, 1997.

3. Wakisaka Y, Okuzono Y, Taenaka Y, et al: Development of a flow estimation and control system of an implantable centrifugal blood pump for circulatory assist. Artif Organs 22: 488–492, 1998.

4. Wakisaka Y, Okuzono Y, Taenaka Y, Chikanari K, Masuzawa T, Takano H: Establishment of flow estimation for an implantable centrifugal blood pump. ASAIO J 43: 659–662, 1997.

5. Wakisaka Y, Okuzono Y, Taenaka Y, et al: Noninvasive pump flow estimation of a centrifugal blood pump. Artif Organs 21: 651–654, 1997.

6. Goodwin GC, Sin KS: Adaptive Filtering Prediction and Control. Englewood Cliffs NJ, Prentice-Hall, 1984.

7. Brooks DE, Goodwin JW, Seaman GVF: Interactions among erythrocytes under shear. J Appl Physiol 28: 172–177, 1970.

8. Kisanuki M, Yoshizawa M, Abe K, et al: Parameter estimation of cardiovascular dynamics for artificial heart control. Jpn J Artif Organs 24: 1099–1106, 1995.

9. Yoshizawa M, Kisanuki M, Abe K, Yambe T, Nitta S: Hemodynamic waveform estimation of artificial hearts using artificial neural networks. Biomed Eng 8: 76–85, 1996.

Copyright © 2002 by the American Society for Artificial Internal Organs