Characterization of a pediatric rotary blood pump

Abstract

Introduction

A ventricular assist device (VAD) is an electromechanical pump used to treat heart failures. For designing the physiological control system for a VAD, one needs a mathematical model and its related parameters. This paper presents a characterization procedure for determining the model parameter values of the electrical, mechanical, and hydraulic subsystems of a pediatric Rotary Blood Pump (pRBP).

Methods

An in vitro test setup consisting of a pRBP prototype, a motor driver module, an acrylic reservoir, mechanical resistance and tubings, pressure and fluid flow sensors, and data acquisition, processing, and visualization system. The proposed procedure requires a set of experimental tests, and a parameter estimation algorithm for determining the model parameters values.

Results

The operating limits of the pRBP were identified from the steady-state data. The relationship between the pressure head, flow rate, and the rotational speed of the pRBP was found from the static tests. For the electrical and mechanical subsystems, the dc motor model has a viscous friction coefficient that varies nonlinearly with the flow. For the hydraulic subsystem, the pressure head is assumed to be a sum of terms related to the resistance, the inertance, the friction coefficient, and the pump speed.

Conclusion

The proposed methodology was successfully applied to the characterization of the pRBP. The combined use of static and dynamic tests provided a precise lumped parameter model for representing the pRBP dynamics. The agreement, regarding mean squared deviation, between experimental and simulated results demonstrates the correctness and feasibility of the characterization procedure.

Keywords
Rotary blood pump; Ventricular assist device; Centrifugal flow pump; Lumped parameter model; System identification

Introduction

Ventricular assist devices are an important and increasingly prevalent treatment for adult and pediatric patients exhibiting advanced heart failure. These devices can be used as a bridge to heart transplant or destination therapy, and they are classified according to the pumping principle as volume-displacement (pulsatile-) or continuous-flow pumps (RBPs) ( Stanfield and Selzman, 2013 ).

The evolution of VADs from pulsatile pumps to rotary blood pumps has started a new era for the treatment of heart failures. The RBP provides a continuous flow of blood for the patient and the increase of its clinical usage comes from their smaller size, improved durability and increased survival with less morbidity when compared to the performance of pulsatile pumps ( Moazami et al., 2013 ).

Typical impeller types for RBPs, such as axial or centrifugal pumps, are designed to operate around a specific speed, pressure, and flow. Classically, axial-flow pumps tend to be used for high flow rates and low-pressure head conditions, whereas centrifugal-flow ones are designed for low flow rates and high-pressure head conditions ( Stanfield et al., 2012 ).

For the development of the physiological control system, the steady-state and dynamic operation of RBPs must be investigated. Within the model-based design framework, rotary blood pumps used as VADs are usually represented by lumped parameter models, named as 0D models ( Capoccia, 2015 ; Moscato et al., 2009 ). This type of model is the most adopted because it uses a relatively small number of parameters to simulate the cardiovascular hemodynamics. Besides, the concept of analog circuitry is used, in which the electric current represents the fluid flow in the blood vessels, and the electrical voltage represents the pressures in the circulatory system ( Shi et al., 2011 ).

The objectives of this paper are to obtain the lumped parameter model and to present a characterization procedure for a pediatric rotary blood pump (pRBP) prototype.

Methods

This section describes the in vitro experimental test setup configuration, the dynamic equations used for representing the behavior of the electrical, mechanical and hydraulic subsystems, and the proposed characterization procedure for pRBP.

In vitro experimental test setup

The in vitro experimental test setup used to characterize the pRBP is shown in Figure 1 . This test setup consists of: a rotary blood pump prototype, a motor driver module, pressure and fluid flow sensors, data acquisition boards, an acrylic reservoir, a mechanical resistance, and tubings. A computer (not shown) running the Microsoft Windows operating system provides the data processing and data visualization. A USB interface provides the data transmission from the test setup to the computer. The structure of the experimental test setup is relatively simple to build and allows to put the pRBP to operate at different pressures and fluid flow levels.

Figure 1
Photograph of the in vitro experimental test setup. The parts identification and description are provided in the white box.

Figure 2 shows a simplified block diagram to explain the coupling of the electrical, mechanical and hydraulic subsystems. The dc motor is the actuator of the electrical subsystem, in which the motor driver module provides the armature voltage. The pump is connected to the reservoir, the pressures and flow may be measured at the inlet and outlet of the pump. It is worth to point out that this type of test setup structure is widely used in RBPs characterization ( Ganushchak et al., 2006 ; Jahanmir et al., 2008 ; Pennings et al., 2013 ).

Figure 2
Block diagram for the in vitro experimental test setup.

In the hydraulic subsystem, the pRBP (20 mL ejection volume) is a centrifugal rotary blood pump to be used as a bedside pump for the patient. In the mechanical subsystem, there is a coupling between the pRBP and the dc motor shaft. The two parts of the coupler, one at dc motor side and the other at the impeller side have four permanent magnets.

The coupling between the motor shaft and the RBP impeller is magnetic since it causes less friction with blood; thus heat generation is minimized, and therefore, less formation of thrombi and occurrence of hemolysis ( Pérez, 2009 ). In this type of pump, the blood flow is generated by the impeller rotational speed. The impeller is enclosed in a polycarbonate housing with central access inlet connector and peripheral access outlet connector.

The rotary pump inlet and outlet accesses are connected to an acrylic reservoir using tubings. The liquid used in the tests is a solution composed of $2/3$ saline solution and $1/3$ glycerin designed to reproduce the blood viscosity. Additionally, there is a mechanical resistance which changes the cross-section of the tubing to simulate obstructions. The test setup is equipped with sensors that allow the measurement of the pressure head $\Delta p\left(t\right)$ (the pressure difference between the pump outlet ( $p_{out}\left(t\right)$ ) and pump inlet ( $p_{in}\left(t\right)$ )), the fluid flow at the outlet tubing $q\left(t\right)$ , the motor armature voltage $v_a\left(t\right)$ , the armature current $i_a\left(t\right)$ , and the mechanical shaft speed $\omega \left(t\right)$ . The pressure signal is measured with a TruWave disposable pressure transducer (Model number PXMK2051), which is based on the principle of strain-gauge and operates in the pressure range from -50 to +300 mmHg ( Edwards Lifesciences Corporation, 2018 ). The fluid flow is measured with a em-tec Digiflow-EXT1 non-invasive flow sensor (Model number 2017103945), using ultrasound scheme with Clamp-On Transducer CT 1/4" x 1/16” and operates in the flow range ±32 L/min ( Germany em-tec GmbH, 2018 ). The rotational speed is measured by using a 15-hole encoder installed at the dc motor side of the magnetic coupler. The analog signals provided by these sensors are acquired by using National Instruments boards (NI USB-6210 and NI myDAQ) and delivered to the computer via the USB interface. These boards allow both single and multiple 16-bit analog-digital conversions. To avoid data loss during the experimental tests a large first-in-first-out buffer is provided. The setup also has software for selecting the structure of the mechanical and hydrodynamic dynamics, and for estimating all the model parameters; the software runs on a desktop computer not shown in Figure 2 .

Pump modeling

The lumped parameter model for the pRBP prototype system is described by a set of three ordinary differential equations. The equivalent circuit representation for this lumped parameter model is as shown in Figure 3 .

Figure 3
Lumped parameter representation of the pRBP.

The electrical subsystem is composed of a dc motor. The differential equation for this subsystem is found by applying the Kirchhoff’s voltage law to the armature circuit. Thus, the model for the electrical subsystem is given by

$\frac{d{i}_a\left(t\right)}{dt}=\frac{1}{L_a}{v}_a\left(t\right)-\frac{R_a}{L_a}{i}_a\left(t\right)-\frac{k_e}{L_a}\omega \left(t\right),$ (1)

where $R_a$ is the armature resistance, $L_a$ is the armature inductance, $k_e$ is the back-emf constant, $e_a\left(t\right)$ is back-emf voltage, and $\omega \left(t\right)$ is the rotor speed, in which $e_a\left(t\right)=k_e\omega \left(t\right)$ .

For the mechanical dynamics, the difference of the electromagnetic torque and pump torque determines the energy used for accelerating or decelerating the impeller. There is an axial-type magnetic coupler between the dc motor and the pump impeller. Permanent magnets are disposed on a disk fixed at the dc motor shaft spinning with rotational speed ${\omega }_m$ . At the pump impeller side, the magnets are disposed over a polycarbonate disk rotating with rotational speed ${\omega }_p$ ( Lubin et al., 2014 ). In this paper, for simplicity, we assume that the coupler is rigid, i.e., ${\omega }_m\left(t\right)={\omega }_p\left(t\right)=\omega \left(t\right)$ , and thus the motion equation is given by

$J_m\frac{d\omega \left(t\right)}{dt}-c\omega \left(t\right)=T_e-T_l,$ (2)

where $J_m$ is the moment of inertia, $c$ is the viscous friction coefficient, $T_e$ is the electromagnetic torque, given by $T_e=k_e{i}_a\left(t\right)$ , and $T_l$ is the load torque (pump torque). Defining $f\left(q\left(t\right),\omega \left(t\right)\right)=c\omega \left(t\right)+T_l$ as the net pump torque one may rewrite the motion equations as

$\frac{d\omega \left(t\right)}{dt}=\frac{k_e}{J_m}{i}_a\left(t\right)-\frac{1}{J_m}f\left(q\left(t\right),\omega \left(t\right)\right)$ . (3)

The reason for writing the motion equation like this is related with the fact that the expression for $f\left(q\left(t\right),\omega \left(t\right)\right)$ is not known a priori, and must be defined in the characterization procedure. The pump torque acting on the dc motor shaft depends on the type of coupling device and on how the impeller blades interact with the flow. In this paper, we consider that $f\left(q\left(t\right),\omega \left(t\right)\right)$ is one out of the five possible structures given below

$\text{Model I}:f\left(q\left(t\right),\omega \left(t\right)\right)=c\omega \left(t\right),$ (4)

$\text{Model II}:f\left(q\left(t\right),\omega \left(t\right)\right)=\left[a{q}^2\left(t\right)+c\right]\omega \left(t\right),$ (5)

$\text{Model III}:f\left(q\left(t\right),\omega \left(t\right)\right)=\left[a{q}^2\left(t\right)+c\right]\omega \left(t\right)+bq\left(t\right){\omega }^2\left(t\right),$ (6)

$\text{Model IV}:f\left(q\left(t\right),\omega \left(t\right)\right)=\left[a{q}^2\left(t\right)+c\right]\omega \left(t\right)+d{\omega }^3\left(t\right),$ (7)

$\text{Model V}:f\left(q\left(t\right),\omega \left(t\right)\right)=\left[a{q}^2\left(t\right)+c\right]\omega \left(t\right)+bq\left(t\right){\omega }^2\left(t\right)+d{\omega }^3\left(t\right)$ . (8)

In Model I, the function $f\left(q\left(t\right),\omega \left(t\right)\right)$ is dependent only on the linear component of the viscous friction ( $c$ ) generated by the dc motor. In Model II, a non-linear term quadratic (a ) with the flow is added to the constant viscous friction coefficient. In Model III, a bi-linear term (b) involving the flow multiplied by the square of the rotational speed is included. In Model IV, the bi-linear term is replaced by one that depends on the cube of the rotational speed (d). In Model V, all the previous terms are combined ( Choi et al., 1997 ; Lim et al., 2007 ).

For the hydraulic dynamics, the pressure head, fluid flow, and rotational speed relationship can be described by using Euler’s pump and turbine equations. However, there is no simple theoretical method for finding this relationship without using finite element analysis. Thus, similarly to what was said for the motion equation, the function that relates the pressure head with fluid flow and rotational speed is not known a priori and also must be defined in the characterization procedure. In this paper, we consider that $\Delta p\left(t\right)$ is one of the five possible models given below

$$\text{Model I}:\Delta p\left(t\right)=R_{p}q\left(t\right)+\beta {\omega }^2\left(t\right),$$ (9)

$\text{Model II}:\Delta p\left(t\right)=R_{p}q\left(t\right)+F_p{q}^2\left(t\right)+\beta {\omega }^2\left(t\right),$ (10)

$\text{Model III}:\Delta p\left(t\right)=R_{p}q\left(t\right)+F_p{q}^2\left(t\right)+\beta {\omega }^2\left(t\right)+J_p\frac{d\omega \left(t\right)}{dt},$ (11)

$\text{Model IV}:\Delta p\left(t\right)=R_{p}q\left(t\right)+F_p{q}^2\left(t\right)+\beta {\omega }^2\left(t\right)+L_p\frac{dq\left(t\right)}{dt},$ (12)

$\text{Model V}:\Delta p\left(t\right)=R_{p}q\left(t\right)+F_p{q}^2\left(t\right)+\beta {\omega }^2\left(t\right)+L_p\frac{dq\left(t\right)}{dt}+J_p\frac{d\omega \left(t\right)}{dt}.$ (13)

In Model I, the pressure head is assumed to be a sum of a term proportional to the flow (related to the pump resistance $R_p$ ), and a nonlinear term quadratic with the rotational speed ( $\beta$ ), in which $e_p\left(t\right)=\beta {\omega }^2\left(t\right)$ . In Model II, non-linear term quadratic ( $F_p$ ) with the flow is added. In Model III, a term proportional to the time derivative of the rotational speed (related to the impeller pump’s moment of inertia $J_p$ ) is included. In Model IV, the term proportional to the time derivative of the rotational speed is replaced by one proportional to the time derivative of the flow (related to the pump’s inertance $L_p$ ). In Model V, all the previous terms are combined ( Pirbodaghi et al., 2011 ; Pirbodaghi, 2017; Zhang et al., 2010 ).

As mentioned before, the selection of the structure for the mechanical and hydrodynamic dynamics, and the estimation of all the model parameters is provided by a software. Thus, the continuous-time models used for representing the electrical, mechanical and hydraulic subsystems must be converted into discrete difference equations, suitable for numerical computing. Among the several methods for deriving the discrete-time representation ( Santina and Stubberud, 2010 ), in this paper, for simplicity, the forward Euler method was used, i.e., $\frac{dx\left(t\right)}{dt}\approx \frac{x\left(t+\Delta t\right)-x\left(t\right)}{\Delta t}$ , where $x\left(t\right)$ may represent $i_a\left(t\right)$ , $\omega \left(t\right)$ , $q\left(t\right)$ . Applying forward Euler method to the electrical, mechanical and hydraulic subsystems models the following discrete-time representations are obtained

$i_a\left(k+1\right)=\left(1-\frac{h{R}_a}{L_a}\right)i_a\left(k\right)+\frac{h}{L_a}{v}_a\left(k\right)-\frac{2h{k}_e}{L_a}\omega \left(k\right),$ (14)

$\omega \left(k+1\right)=\left(1-\frac{hc}{J_m}\right)\omega \left(k\right)+\frac{h{k}_e}{J_m}{i}_a\left(k\right)-\frac{ha}{J_m}{q}^2\left(k\right)\omega \left(k\right)-\frac{hb}{J_m}q\left(k\right){\omega }^2\left(k\right)-\frac{hd}{J_m}{\omega }^3\left(k\right),$ (15)

$q\left(k+1\right)=\left(1-\frac{h{R}_p}{L_p}\right)q\left(k\right)-\frac{h{F}_p}{L_p}{q}^2\left(k\right)-\frac{h\beta }{L_p}{\omega }^2\left(k\right)+\frac{h}{L_p}\Delta p\left(k\right)-\frac{J_p}{L_p}\left[\omega \left(k+1\right)-\omega \left(k\right)\right],$ (16)

where $k=\mathrm{1,}\cdots ,N$ , with $N$ is the number of data samples and $h=\Delta t$ is the sampling interval. In the case of the mechanical and hydraulic models, we only show the equations for the fifth cases (Model V), since from these expressions one may derive the other representations (Models I-IV) just by deleting the terms not being used.

Pump characterization procedure

The selection of the structure for the mechanical and hydrodynamic dynamics and the estimation of all the model parameters of the pRBP is based on a sequence of experimental tests and a numerical algorithm. The proposed characterization procedure for this purpose is sketched in the flowchart of Figure 4 . In the following a detailed explanation of all steps of the characterization procedure will be provided.

Figure 4
Flow chart for description of the pump characterization procedure. In this case, each interaction is represented in Roman numerals to indicate the verification of each proposed model.

Definition of operating limits

An initial test using the pRBP is executed for defining the operating limits, regarding the maximum and minimum values of the signals available on the in vitro experimental setup. This step is performed as follows:

1. a

   With the mechanical resistance fully open (see tag 2 at Figure 1 ), the armature voltage (see $v_a\left(t\right)$ at Figure 2 ) is changed from 0 V to the maximum voltage provided by motor driver module in steps of 1 V at a time; at each step the voltage level is kept constant for 1 min;

2. b

   After 1 min, when the steady-state is reached, the armature voltage, the armature current, the rotational speed, the fluid flow and the pressures (at the pump inlet and pump outlet) are recorded;

3. c

   Repeat (a) and (b) until the magnetic coupler decouples the motor shaft from the pump. The decoupling is detected by a sudden reduction of the armature current, rotational speed, and fluid flow.

Static test

The traditional method for characterization of rotary pumps consists in generating pressure head x flow plots, i.e., $\Delta p\times q$ curves ( Ganushchak et al., 2006 ). In this method, the flow rate, and the pressure head at different rotational speeds, under steady-state conditions, are recorded. The procedure to obtain hydrodynamic performance is as follows:

1. i

   In the motor driver module, the pump speed is changed from minimum speed (the value that guarantees the complete circulation of fluid in the experimental setup) to the maximum speed (the value that allows pumping without to sink the flow during the tubing obstruction) by applying a 500 RPM step.

2. ii

   For each rotational speed, the flow is changed in steps of 1 L/min and maintaining each step for 1 min based on the mechanical resistance, from fully closed to fully open as shown in Figure 5 a.

   

   Figure 5
   Representation of the (a) $q\left(t\right)$ and (b) $\Delta p$ (t) steps in the static test; (c) $\Delta p\times q$ resulting curves for $n$ different rotational speeds.
   

3. iii

   After stabilizing the flow at each step, the inlet and outlet pump pressures are recorded as shown in Figure 5 b. At the end of the test, $\Delta p\times q$ curves are plotted with the measured values, as shown in Figure 5 c.

Dynamic test

To select the structure for the mechanical and hydrodynamic dynamics and to estimate all the pRBP model parameters, a dynamic test is done on the in vitro experimental test setup based on operating limits already found. The procedure used for this test is as follows:

1. i

   With the mechanical resistance fully open, a variable frequency voltage signal is applied to the armature of dc motor for 1 min.

2. ii

   During this 1 min, the armature voltage, the armature current, the rotational speed, the flow and the inlet/outlet pressures are recorded. At the end of the test, the recorded data is split into two parts, one for parameter estimation and another part for model validation.

Parameter estimation

Based on the discrete-time models given by Equations 14 , 15 and 16 one may derive one-step-ahead predictors for formulating the parameter estimation problems ( Ljung, 1987 ). For this purpose the discrete-time models must be rewritten as

$\widehat{y}(k|{\theta }_i)={\phi }^T\left(k\right){\theta }_i$ . (17)

Based on $Z_i^{N_E}$ , the experimental data matrix recorded during the dynamic test, a cost function $V\left({\theta }_i,Z_i^{N_E}\right)$ is defined. The cost function is expressed by

$$V\left({\theta }_i,Z_i^{N_E}\right)=\frac{1}{N_E}{\displaystyle \sum _{k=1}^{N_E}{\left[y_i\left(k\right)-\widehat{y}\left(k\text{|}{\theta }_i\right)\right]}^2},i\in \left\{\text{e},\text{m},\text{h}\right\},$$ (18)

where $N_E$ is the number of data points selected for parameter estimation, $i$ denotes the electrical (e), mechanical (m) or hydraulic (h) discrete-time subsystem model. The parameter estimate for each subsystem model is found with a numerical algorithm designed to solve ( Ljung, 1987 )

${\widehat{\theta }}_i=\underset{{\theta }_i}{\text{argmin}}V\left({\theta }_i,Z_i^{N_E}\right),i\in \left\{\text{e},\text{m},\text{h}\right\}$ . (19)

Model validation

For determining the quality of the estimated model, the part of the data matrix not used in Equation 19 is used to compute a figure of merit, the root mean square error (RMSE) index. This figure of merit indicates the accuracy of the estimated model.

$W{}_d\left({\widehat{\theta }}_i^j,Z_i^{N_V}\right)=\sqrt{\frac{1}{N_V}{\displaystyle \sum _{k=1}^{N_V}{\left[y_i\left(k\right)-\widehat{y}\left(k\text{|}{\widehat{\theta }}_i^j\right)\right]}^2}},i\in \left\{\text{e},\text{m},\text{h}\right\},j\in \left\{\text{I},\text{II},\text{III},\text{IV},\text{V}\right\},$ (20)

where $N_V=N-N_E$ .

The block named “Model choice” (see Figure 4 ) is implemented by using the above figure of merit. Once all the parameters of all models (Models I-V for the mechanical and for the hydraulic subsystems) are available, one computes ${\mathcal{W}}_d\left({\widehat{\theta }}_m\right)=\left\{W{}_d\left({\widehat{\theta }}_m^I\right),W{}_d\left({\widehat{\theta }}_m^{II}\right),W{}_d\left({\widehat{\theta }}_m^{III}\right),W{}_d\left({\widehat{\theta }}_m^{IV}\right),W{}_d\left({\widehat{\theta }}_m^V\right)\right\}$ and ${\mathcal{W}}_d\left({\widehat{\theta }}_h\right)=\left\{W{}_d\left({\widehat{\theta }}_h^I\right),W{}_d\left({\widehat{\theta }}_h^{II}\right),W{}_d\left({\widehat{\theta }}_h^{III}\right),W{}_d\left({\widehat{\theta }}_h^{IV}\right),W{}_d\left({\widehat{\theta }}_h^V\right)\right\}$ and look for the model exhibiting the lowest $W{}_d\left({\widehat{\theta }}_i^j\right),i\in \left\{\text{m},\text{h}\right\},j\in \left\{\text{I},\text{II},\text{III},\text{IV},\text{V}\right\}$ .

Rebuilding of static curves

An additional step, usually not considered in parameter estimation problems is introduced in the proposed characterization procedure. This additional step consists in using the model issued form the parameter estimation algorithm to simulate the static test and compare the simulated results with the actual data recorded previously. The final quality of the model is defined regarding another figure of merit. This figure of merit is expressed by

$$W{}_{s,\omega }\left({\widehat{\theta }}_h^j\right)=\sqrt{\frac{1}{M}{\displaystyle \sum _{k=1}^M{\left[q_s\left(k\right)-{\widehat{q}}_s\left(k\text{|}{\widehat{\theta }}_h^j\right)\right]}^2}},\omega \in \left[{\omega }_1,{\omega }_n\right],j\in \left\{\text{I},\text{II},\text{III},\text{IV},\text{V}\right\},$$ (21)

where $M$ is the number of data points recorded in the static test, $q_s\left(k\right)$ is the actual fluid flow measurement, and ${\widehat{q}}_s(k|{\widehat{\theta }}_{h_i})$ is the simulated fluid flow computed with estimated parameter set at the same operating conditions of the static test. This figure of merit is like a cross-validation procedure that corroborates the correctness of the proposed approach specifically regarding the hydrodynamic behavior.

Results

In this section, the data collected at the three first steps of the proposed characterization procedure are presented, compared and discussed.

Operating limits

From the measurements collected at the in vitro experimental setup, the operating limits of the pRBP were defined. In this case, the voltage applied to the armature of the dc motor can keep the pump in full operation up to 21 V, with the current saturating at 2 A and the rotational speed around 6000 RPM, as can be seen in Figure 6 a. After this operating point, the flow, and the inlet/outlet pressures goes to zero; this fact is an indication that there is an uncoupling between dc motor and the pump impeller, as can be seen in Figure 6 b.

Figure 6
Definition of the pump operating limits for: (a) electrical and mechanical quantities; (b) hydraulic quantities.

Static test

After defining the operating range, the static test was done using the pRBP. In this work, the minimum rotational speed for start pumping is 1000 RPM, whereas the maximum rotational speed, with the mechanical resistance fully closed, is 4000 RPM; above that value, the pump collapses. Thus, in all the subsequent tests the pump speed was restricted to the interval [1000,4000] RPM. Figure 7 shows the plots of $q\left(t\right)\times t$ and $\Delta p\left(t\right)\times t$ as obtained in one of the static tests. From the signals recorded in all the static tests, the performance curves were plotted as shown in Figure 8 .

Figure 7
The $q\left(t\right)$ and $\Delta p\left(t\right)$ steps used in the static test for $\omega =3500$ RPM.

Figure 8
The $\Delta p\times q$ measured curves for the pRBP.

Dynamic test

The dynamic test provides the measurements of armature voltage, stator current, flow, pressure head and rotational speed available on the experimental test setup to be used in the parameter estimation and validation tasks. As mentioned before the experimentally collect data are split into two parts. Eighty percent of data were used for parameter estimation, whereas twenty percent of them were used for model validation. In both cases, the units of measurement of flow and rotational speed were converted to standard units, i.e., ml/s and rad/s, respectively.

For estimating the discrete-time subsystems model parameters, the extended recursive least squares (ERLS) algorithm ( Åström and Wittenmark, 1994 ; Chen, 2004 ) was used; it allows one to estimate numerical values for $R_a$ , $L_a$ , $k_e$ , $J_m$ , $a$ , $b$ , $c$ , $d$ , $R_p$ , $F_p$ , $\beta$ , $L_p$ and $J_p$ . In this algorithm, the model is considered as a stochastic system, wherein each differential equation is rewritten as a autoregressive-moving average with exogenous terms (ARMAX) process ( Ljung, 1987 ), according in Equation (22) .

$A\left(q\right)y\left(k\right)=B\left(q\right)u\left(k\right)+C\left(q\right)e\left(k\right),$ (22)

where $y\left(k\right)$ is the model output, $u\left(k\right)$ is the model input, $e\left(k\right)$ is the system noise, $A\left(q\right)=1+a_1{q}^{-1}+\mathrm{...}+a_{n_a}{q}^{-n_a}$ , $B\left(q\right)=b_1{q}^{-1}+\mathrm{...}+b_{n_b}{q}^{-n_b}$ and $C\left(q\right)=1+c_1{q}^{-1}+\mathrm{...}+c_{n_c}{q}^{-n_c}$ are the polynomials in the forward shift operator $q$ , with orders $n_a$ , $n_b$ and $n_c$ , respectively. The ARMAX model cannot be converted directly to a regression model, since the noise variables $e\left(t\right)$ are not known a priori. Thus, to derive a regression model suitable approximation, in which the noise is approximated by the residual vector $\epsilon \left(k\right)$ , as expressed in Equation (23) .

$\epsilon \left(k\right)=y\left(k\right)-{\phi }^T\left(k\right)\widehat{\theta }\left(k\right),$ (23)

where $\theta ={[a_1\mathrm{...}{a}_{n_a}{b}_1\mathrm{...}{b}_{n_b}{c}_1\mathrm{...}{c}_{n_c}]}^T$ is the parameters vector and ${\phi }^T\left(k\right)=\left[-y\left(k-1\right)\mathrm{...}-y\left(k-n\right)u\left(k-1\right)\mathrm{...}u\left(k-n\right)\epsilon \left(k-1\right)\mathrm{...}\epsilon \left(k-n\right)\right]$ is the regression vector. For the discrete-time subsystems model, the parameters vectors according to the ARMAX structure were defined as given in Equations 24 - 34 , with $h=0.01$ seconds.

• Electrical dynamics:

  ${\theta }_e={\left[\left(1-\frac{h{R}_a}{L_a}\right)\frac{h}{L_a}\frac{2h{k}_e}{L_a}\right]}^T$(24)

• Mechanical dynamics:

  ${\theta }_m^I={\left[\left(1-\frac{hc}{J_m}\right)\frac{h{k}_e}{J_m}\right]}^T$(25)
  ${\theta }_m^{II}={\left[\left(1-\frac{hc}{J_m}\right)\frac{h{k}_e}{J_m}\frac{ha}{J_m}\right]}^T$(26)
  ${\theta }_m^{III}={\left[\left(1-\frac{hc}{J_m}\right)\frac{h{k}_e}{J_m}\frac{ha}{J_m}\frac{hb}{J_m}\right]}^T$(27)
  ${\theta }_m^{IV}={\left[\left(1-\frac{hc}{J_m}\right)\frac{h{k}_e}{J_m}\frac{ha}{J_m}\frac{hd}{J_m}\right]}^T$(28)
  ${\theta }_m^V={\left[\left(1-\frac{hc}{J_m}\right)\frac{h{k}_e}{J_m}\frac{ha}{J_m}\frac{hb}{J_m}\frac{hd}{J_m}\right]}^T$(29)

• Hydraulic dynamics:

  ${\theta }_h^I={\left[\frac{\beta }{R_p}\frac{1}{R_p}\right]}^T$(30)
  ${\theta }_h^{II}={\left[\frac{F_p}{R_p}\frac{\beta }{R_p}\frac{1}{R_p}\right]}^T$(31)
  ${\theta }_h^{III}={\left[\frac{F_p}{R_p}\frac{\beta }{R_p}\frac{1}{R_p}\frac{J_p}{R_p}\right]}^T$(32)
  ${\theta }_h^{IV}={\left[\left(1-\frac{h{R}_p}{L_p}\right)\frac{F_p}{L_p}\frac{\beta }{L_p}\frac{1}{L_p}\right]}^T$(33)
  ${\theta }_h^V={\left[\left(1-\frac{h{R}_p}{L_p}\right)\frac{F_p}{L_p}\frac{\beta }{L_p}\frac{1}{L_p}\frac{J_p}{L_p}\right]}^T$(34)

These vectors have been used in the ERLS algorithm. At the of the estimation step, by using the estimated parameters, the cost functions, and the RMSE indexes have been computed. The results are presented in Tables 1 , 2 and 3 .

Table 1
Estimated parameters for 0D model - electrical subsystem.

${\widehat{R}}_a$	${\widehat{L}}_a$	${\widehat{k}}_e$	$V\left({\theta }_e\right)$	$W{}_d\left({\widehat{\theta }}_e\right)$
( $\Omega$ )	( $mH$ )	$\left(\frac{V\cdot s}{rad}\right)$	(A)	(A)
1.87	38.82	$1.24\times {10}^{-2}$	$4.06\times {10}^{-3}$	$6.24\times {10}^{-5}$

Table 2
Estimated parameters for 0D model - mechanical subsystem.

Model	${\widehat{J}}_m$	$â$	$\widehat{b}$	$\widehat{c}$	$\widehat{d}$	$V\left({\theta }_m\right)$	$W{}_d\left({\widehat{\theta }}_m\right)$
	$\left(kg\cdot m{}^2\right)$	$\left(\frac{kg\cdot m{}^2\cdot s}{m{L}^2\cdot rad}\right)$	$\left(\frac{kg\cdot m{}^2\cdot s}{mL\cdot rad{}^2}\right)$	$\left(\frac{kg\cdot m{}^2}{rad\cdot s}\right)$	$\left(\frac{kg\cdot m{}^2\cdot s}{rad{}^3}\right)$	$\left(\frac{rad}{s}\right)$	$\left(\frac{rad}{s}\right)$
I	$1.03\times {10}^{-4}$	-	-	$3.92\times {10}^{-5}$	-	0.28	$4.00\times {10}^{-4}$
$\text{*}$ II	$9.38\times {10}^{-5}$	- $2.90\times {10}^{-9}$	-	$6.76\times {10}^{-5}$	-	0.28	$3.97\times {10}^{-4}$
III	$1.03\times {10}^{-4}$	- $2.47\times {10}^{-8}$	$6.56\times {10}^{-9}$	$1.58\times {10}^{-5}$	-	0.33	$4.49\times {10}^{-4}$
IV	$1.03\times {10}^{-4}$	- $1.13\times {10}^{-8}$	-	$1.50\times {10}^{-5}$	$8.04\times {10}^{-10}$	0.33	$4.51\times {10}^{-4}$
V	$1.03\times {10}^{-4}$	$3.62\times {10}^{-7}$	- $1.80\times {10}^{-7}$	$1.48\times {10}^{-5}$	$2.26\times {10}^{-8}$	0.34	$4.59\times {10}^{-4}$

Table 3
Estimated parameters for 0D model - hydraulic subsystem.

Model	${\widehat{R}}_p$	${\widehat{F}}_p$	$\widehat{\beta }$	${\widehat{L}}_p$	${\widehat{J}}_p$	$V\left({\theta }_h\right)$	$W{}_d\left({\widehat{\theta }}_h\right)$
	$\left(\frac{mmHg\cdot s}{mL}\right)$	$\left(\frac{mmHg\cdot s{}^2}{mL{}^2}\right)$	$\left(\frac{mmHg\cdot s{}^2}{rad{}^2}\right)$	$\left(\frac{mmHg\cdot s{}^2}{mL}\right)$	$\left(\frac{mmHg\cdot s{}^2}{rad}\right)$	$\left(\frac{mL}{s}\right)$	$\left(\frac{mL}{s}\right)$
I	7.88	-	- $3.71\times {10}^{-3}$	-	-	5.27	$3.53\times {10}^{-3}$
II	8.06	-0.04	- $9.51\times {10}^{-4}$	-	-	1.06	$2.23\times {10}^{-3}$
III	7.43	-0.04	- $9.89\times {10}^{-4}$	-	0.07	0.95	$1.91\times {10}^{-3}$
$\text{*}$ IV	2.51	0.01	- $1.14\times {10}^{-3}$	0.68	-	1.04	$1.92\times {10}^{-4}$
V	2.47	0.01	- $1.11\times {10}^{-3}$	0.69	- $5.16\times {10}^{-5}$	1.05	$1.93\times {10}^{-4}$

In Table 1 , the parameters calculated for electrical subsystem are suitable for the model because these parameters yielded a small value for $V\left({\theta }_e\right)$ and $W{}_d\left({\widehat{\theta }}_e\right)$ . Although the Models I and II have presented the same $V\left({\theta }_m\right)$ in Table 2 , Model II was chosen for mechanical dynamics because it exhibted the lowest $W{}_d\left({\widehat{\theta }}_m\right)$ .

As seen in Table 3 , the Model III was chosen as the initial model for hydraulic dynamics with minimum $V\left({\theta }_h\right)$ , but the estimated static curves did not track the reference curves, represented by the high value of the $W{}_d\left({\widehat{\theta }}_h\right)$ . In sequence, the Model IV was analyzed, in which the estimated curves provided a good representation of the experimental data with lowest $W{}_d\left({\widehat{\theta }}_h\right)$ .

Furthermore, the comparison between the measured and estimated signals for the electrical model, mechanical Model II and hydraulic Model IV, are shown in Figures 9 99c, respectively.

Figure 9
(a) Comparing $i_a\left(t\right)$ and ${\widehat{i}}_a\left(t\right)$ ; (b) Comparing $\omega \left(t\right)$ and $\widehat{\omega }\left(t\right)$ ; (c) Comparing $q\left(t\right)$ and $\widehat{q}\left(t\right)$ .

Thus, the static performance curves obtained experimentally and its rebuilding using the hydraulic Model IV can be observed in Figure 10 .

Figure 10
The $\Delta p\times q$ and $\Delta p\times \widehat{q}$ curves for the pRBP.

Discussion

With the experimental test setup and the characterization procedure, it was possible to find a suitable model for pRBP prototype. The pRBP dynamic model is composed of three coupled differential equations, the first one for the electrical subsystem, the second for the mechanical subsystem and the third one for the hydraulic subsystem. For the parameter estimation, the ERLS algorithm was adopted. The electrical parameters are given in Table 1 , the mechanical parameters are given in Table 2 (line marked with a star), and the hydraulic parameters are given in Table 3 (line marked with a star).

Besides that, the pump characterization procedure provided a simple methodology for determining the values of the continuous-time model parameters requiring relatively few a priori information, specifically in respect to the behavior of the mechanical and hydraulic subsystems. It is important to point out that the proposed characterization procedure is quite generic and can be applied when there are structural changes in the experimental setup, for example, the use a different type of motor or a pump with a different impeller type.

Also, it has been recognized, during the static tests as described in ( Ganushchak et al., 2006 ), that if this pRBP operates in the area to the right of the operating zone (in this work, $\omega >$ 4000 RPM) causes extreme turbulence and cavitation. Even if this pRBP operates in the area to the left of the zone (in this work, $\omega <$ 1000 RPM), its operation is very inefficient, with excessive recirculation of fluid inside the pump and cavitation.

The application of the proposed methodology to a new pump where the motor is a BLDC type is underway. On the other hand, the connection of the pRBP prototype to an in vitro cardiovascular simulator is also underway.

Acknowledgements

The authors thank COPELE-PPgEE, CNPq, FAPESP Grant 201250283-6 and FINEP Grant 01.14.0177.00 for the financial support to the development of this project.

• How to cite this article: Melo TR, Vasconcelos FJS, Ribeiro LHRD, Bacht S, Cestari IA, Rocha Neto JS, Lima AMN. Characterization of a pediatric rotary blood pump. Res Biomed Eng. 2018. DOI: 10.1590//2446-4740.180046

References

• Åström KJ, Wittenmark B. Adaptive control. 2nd ed. Boston: Addison-Wesley Longman Publishing Co., Inc.; 1994.
• Capoccia M. Development and characterization of the arterial windkessel and its role during left ventricular assist device assistance. Artif Organs. 2015; 39(8):E138-53. http://dx.doi.org/10.1111/aor.12532. PMid:26147912.
  » http://dx.doi.org/10.1111/aor.12532
• Chen H. Extended recursive least squares algorithm for nonlinear stochastic systems. In: Proceedings of the American Control Conference; 2004 June 30-July 2; Boston, Massachusetts, USA. USA: IEEE Control Systems Society; 2004. p. 4758-63. http://dx.doi.org/10.23919/ACC.2004.1384065.
  » http://dx.doi.org/10.23919/ACC.2004.1384065
• Choi S, Boston JR, Thomas D, Antaki JF. Modeling and identification of an axial flow blood pump. In: Proceedings of the 1997 American Control Conference (Cat. No.97CH36041); 1997 Jun 6-6; Albuquerque, New Mexico, USA. USA: IEEE Control Systems Society; 1997. p. 3714–5. http://dx.doi.org/10.1109/ACC.1997.609538.
  » http://dx.doi.org/10.1109/ACC.1997.609538
• Edwards Lifesciences Corporation. Truwave disposable pressure transducer [Internet]. Irvine: USA Edwards Lifesciences LLC; 2018. [cited 2018 June 16]. Available from: http://www.edwards.com/eu/products/pressuremonitoring/pages/truwave.aspx
  » http://www.edwards.com/eu/products/pressuremonitoring/pages/truwave.aspx
• Ganushchak Y, van Marken Lichtenbelt W, van der Nagel T, de Jong DS. Hydrodynamic performance and heat generation by centrifugal pumps. Perfusion. 2006; 21(6):373-9. http://dx.doi.org/10.1177/0267659106074003. PMid:17312862.
  » http://dx.doi.org/10.1177/0267659106074003
• Germany em-tec GmbH. Experts in non-invasive flow measurement [Internet] [cited 2018 June 16]. Berlin: Germany em-tec GmbH; 2018. Available from: http://www.em-tec.com/
  » http://www.em-tec.com/
• Jahanmir S, Hunsberger AZ, Heshmat H, Tomaszewski MJ, Walton JF, Weiss WJ, Lukic B, Pae WE, Zapanta CM, Khalapyan TZ. Performance characterization of a rotary centrifugal left ventricular assist device with magnetic suspension. Artif Organs. 2008; 32(5):366-75. http://dx.doi.org/10.1111/j.1525-1594.2008.00559.x. PMid:18471166.
  » http://dx.doi.org/10.1111/j.1525-1594.2008.00559.x
• Lim E, Cloherty SL, Reizes JA, Mason DG, Salamonsen RF, Karantonis DM, et al. A dynamic lumped parameter model of the left ventricular assisted circulation. In: Proceedings of the 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society; 2007 Aug 22-26; Cité Internationale, Lyon, France. USA: IEEE Engineering in Medicine and Biology Society; 2007. p. 3990-3. http://dx.doi.org/10.1109/IEMBS.2007.4353208. PMid: 18002874.
  » http://dx.doi.org/10.1109/IEMBS.2007.4353208
• Ljung L. System identification: theory for the user. 2nd ed. New Jersey: Prentice Hall; 1987.
• Lubin T, Mezani S, Rezzoug A. Experimental and theoretical analysis of axial magnetic coupling under steady-state and transient operation. IEEE Trans Ind Electron. 2014; 61(8):4356-65. http://dx.doi.org/10.1109/TIE.2013.2266087.
  » http://dx.doi.org/10.1109/TIE.2013.2266087
• Moazami N, Fukamachi K, Kobayashi M, Smedira NG, Hoercher KJ, Massiello A, Lee S, Horvath DJ, Starling RC. Axial and centrifugal continuous-flow rotary pumps: A translation from pump mechanics to clinical practice. J Heart Lung Transplant. 2013; 32(1):1-11. http://dx.doi.org/10.1016/j.healun.2012.10.001. PMid:23260699.
  » http://dx.doi.org/10.1016/j.healun.2012.10.001
• Moscato F, Danieli GA, Schima H. Dynamic modeling and identification of an axial flow ventricular assist device. Int J Artif Organs. 2009; 32(6):336-43. http://dx.doi.org/10.1177/039139880903200604. PMid:19670185.
  » http://dx.doi.org/10.1177/039139880903200604
• Pennings KA, Martina JR, Rodermans BF, Lahpor JR, van de Vosse FN, de Mol BA, Rutten MC. Pump flow estimation from pressure head and power uptake for the Heartassist5, Heartmate II, and Heartware VADs. ASAIO J. 2013; 59(4):420-6. http://dx.doi.org/10.1097/MAT.0b013e3182937a3a. PMid:23820282.
  » http://dx.doi.org/10.1097/MAT.0b013e3182937a3a
• Pérez CM. Bombas centrífugas como asistencia ventricular: estado actual. Cir Cardiov. 2009; 16(2):119-24. http://dx.doi.org/10.1016/S1134-0096(09)70156-1.
  » http://dx.doi.org/10.1016/S1134-0096(09)70156-1
• Pirbodaghi T, Weber A, Carrel T, Vandenberghe S. Effect of pulsatility on the mathematical modeling of rotary blood pumps. Artif Organs. 2011; 35(8):825-32. http://dx.doi.org/10.1111/j.1525-1594.2011.01276.x. PMid:21793862.
  » http://dx.doi.org/10.1111/j.1525-1594.2011.01276.x
• Pirbodaghi T. Mathematical modeling of rotary blood pumps in a pulsatile in vitro flow environment. Artif Organs. 2017; 41(8):710-6. http://dx.doi.org/10.1111/aor.12860. PMid:28097669.
  » http://dx.doi.org/10.1111/aor.12860
• Santina M, Stubberud AR. Discrete-time equivalents of continuous-time systems. In: Levine WS, editor. The control handbook: control system fundamentals. Boca Raton: CRC Press; 2010. p. 12.1-12.34.
• Shi Y, Lawford P, Hose R. Review of zero-D and 1-D models of blood flow in the cardiovascular system. Biomed Eng Online. 2011; 10(33):1-38. http://dx.doi.org/10.1186/1475-925X-10-33. PMid:21521508.
  » http://dx.doi.org/10.1186/1475-925X-10-33
• Stanfield JR, Selzman CH, Pardyjak ER, Bamberg S. Flow characteristics of continuous-flow left ventricular assist devices in a novel open-loop system. ASAIO J. 2012; 58(6):590-6. http://dx.doi.org/10.1097/MAT.0b013e31826dcbd9. PMid:22990285.
  » http://dx.doi.org/10.1097/MAT.0b013e31826dcbd9
• Stanfield JR, Selzman CH. In vitro pulsatility analysis of axial-flow and centrifugal-flow left ventricular assist devices. J Biomech Eng. 2013; 135(3):0345051-6. http://dx.doi.org/10.1115/1.4023525. PMid:24231821.
  » http://dx.doi.org/10.1115/1.4023525
• Zhang XT, AlOmari AH, Savkin AV, Ayre PJ, Lim E, Salamonsen RF, et al. In vivo validation of pulsatile flow and differential pressure estimation models in a left ventricular assist device. In: Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology; 2010 31 Aug-4 Sept; Buenos Aires, Argentina. USA: IEEE Engineering in Medicine and Biology Society; 2010. p. 2517-20. http://dx.doi.org/10.1109/IEMBS.2010.5626876.
  » http://dx.doi.org/10.1109/IEMBS.2010.5626876

Publication Dates

• Publication in this collection
  29 Nov 2018
• Date of issue
  Oct-Dec 2018

History

• Received
  16 June 2018
• Accepted
  15 Oct 2018