Numerical Optimal Control of Turbo Dynamic Ventricular Assist Devices

Abstract

The current paper presents a methodology for the derivation of optimal operating strategies for turbo dynamic ventricular assist devices (tVADs). In current clinical practice, tVADs are typically operated at a constant rotational speed, resulting in a blood flow with a low pulsatility. Recent research in the field has aimed at optimizing the interaction between the tVAD and the cardiovascular system by using predefined periodic speed profiles. In the current paper, we avoid the limitation of using predefined profiles by formulating an optimal-control problem based on a mathematical model of the cardiovascular system and the tVAD. The optimal-control problem is solved numerically, leading to cycle-synchronized speed profiles, which are optimal with respect to an arbitrary objective. Here, an adjustable trade-off between the maximization of the flow through the aortic valve and the minimization of the left-ventricular stroke work is chosen. The optimal solutions perform better than constant-speed or sinusoidal-speed profiles for all cases studied. The analysis of optimized solutions provides insight into the optimized interaction between the tVAD and the cardiovascular system. The numerical approach to the optimization of this interaction represents a powerful tool with applications in research related to tVAD control. Furthermore, patient-specific, optimized VAD actuation strategies can potentially be derived from this approach.

1. Introduction

Ventricular assist devices (VADs) are mechanical systems aimed at supporting the blood circulation in patients with severe heart failure. After five decades of intensive development, various devices have become increasingly applied, typically as a “bridge to transplantation”, i.e., for patients with no other alternative left [1]. In addition, by assisting and, thus, unloading the failing cardiac ventricle, VADs have also been shown to contribute to some degree to myocardial recovery [2,3]. Currently, second and third generation VADs are poised to replace the first generation volume displacement VADs. They have been found to be superior in terms of survival rate and reduced adverse events [4].

Feedback control strategies to physiologically adapt the operation of turbo dynamic VADs (tVADs) to the oxygen demand and to optimize left-ventricular unloading have been subject to continuous research [5,6,7,8,9,10]. A number of approaches for realizing a speed modulation of tVADs has been analyzed earlier. A square-wave speed profile was applied by Bearnson et al. [11] and Bourque et al. [12] to increase arterial pulse pressure, where the speed profile was not synchronized to the cardiac cycle. Different types of speed profiles, synchronized to the natural cardiac cycle, have been applied in silico and in vitro to analyze their influence on perfusion, pulse pressure and ventricular unloading [13,14,15,16,17,18]. In vivo experiments have been conducted to validate this approach [19,20]. Common among all previous research on this topic was the preselection of a certain speed profile, such as a sine wave or a square wave. In these studies, the parameters of the selected profile were varied to analyze the effect on relevant physiological quantities. Naturally, the outcome depended on the preselected speed profile.

We present a methodology that circumvents the limitation of preselecting a certain speed profile by the application of model-based numerical optimization. A physiologically motivated objective function was constructed, and the speed profile that minimizes this objective function was found by solving an optimal-control problem (OCP). By defining an OCP, no assumptions about the speed profile needed to be made in advance. Instead, the tVAD is expected to optimally interact with the circulatory system. Due to this absence of restrictions, we hypothesize that the optimized speed trajectory performs better than strategies with either a constant speed or predefined speed profiles.

Optimal control was applied to VADs in previous work. The flow field in a pulsatile volume-displacement VAD was optimized in order to reduce hemolysis and thrombus formation in the device in [21]. In two other studies [22,23], the power consumption of a pusher plate-type VAD was minimized using the optimal control of a linear and time-invariant model of the VAD and a part of the systemic circulation. The methodology described in the current paper differs from previous work as follows: the cardiovascular system, including the beating heart, was entirely integrated in the optimization procedure. Considering this coupled system as a whole entity allowed for the optimization of the interaction between the VAD and the circulation system, rather than of just the performance of the isolated VAD.

In the current paper, a state-space model of the human circulation [24] was used in connection with a mathematical model of a turbo dynamic blood pump. This full model was used for simulations with a constant pump speed, as well as with various sinusoidal speed profiles. A reduced version of the model was used to formulate the OCP, which was solved numerically. The validation of the optimized speed profiles on a hybrid mock circulation [25] indicated that they can potentially be evaluated in a real blood pump.

2. Methods

In this section, three different methods for obtaining a pump speed profile are presented. Either a constant speed, a sinusoidal speed profile synchronized to the cardiac cycle or a speed profile obtained from the numerical solution of an OCP was applied. The first two cases are detailed in Section 2.1, and the formulation, as well as the numerical solution procedure of the OCP are described in Section 2.2. The solutions resulting from these methods were compared using a physiologically motivated objective function in Section 3. For the first two methods, this objective function must be evaluated a posteriori, because the speed profile was predefined and an exhaustive parameter sweep was performed. In the third method, the objective function is a part of the OCP, and its value, as well as the corresponding speed profile result from the solution of this OCP. To ensure a fair comparison of the various pump-actuation strategies, one additional parameter was fixed, namely the total cardiac output. This parameter represents the total amount of blood pumped to the circulatory system by both the heart and the blood pump. The reason for fixing the total cardiac output is outlined in the following.

The autoregulatory mechanisms ensure that the supply of oxygenated blood to the body conforms to its current needs. In order to allow for a fair comparison between any two speed profiles, the total cardiac output has to be equal for both. In the lumped-parameter models used in the current study, there is a unique relation between the activity of the baroreflex mechanisms, defining the arterial resistance, and the mean arterial pressure. These two quantities, namely the arterial resistance and mean pressure, define the total cardiac output. It is therefore sufficient to prescribe the cardiac output to ensure the comparability of different solutions. A value of 5Lmin 5Lmin was imposed, which represents a physiologically plausible value for an adult person at rest [26]. In the model used, this value implied a mean systemic arterial pressure of 100mmHg 100mmHg . The results presented could also be reproduced for other perfusion levels. The fixation of the total cardiac output to 5Lmin 5Lminis also justified by the fact that the current study did not investigate the effect of the VAD actuation on the perfusion, but was focused on the heart-VAD interaction.

The model of the cardiovascular system used here was originally published in [24]. An implementation in MATLAB/Simulink^® (The Mathworks Inc., Natick, MA, USA) of this model was used earlier [25]. The pathologic case defined in [24] was used for all investigations. It is defined by a raised heart rate of 90bpm 90bpm and a reduced ventricular contractility. Figure 1 illustrates the structure of the model. The full model comprises the left and right atria and ventricles, the four heart valves and the systemic and the pulmonary circulation, each with an arterial and a venous section. The four contracting chambers of the heart are modeled based on the principle of time-varying elastance, which is the inverse of the compliance or the stiffness of the heart muscle wall. This stiffness varies during the cardiac cycle, thereby causing a contraction. The time-varying elastance is modeled by specific, empirically-derived curves that are described in the literature [27]. Three independent autoregulatory mechanisms were modeled: the first mechanism adapts the flow resistance in the systemic arterial system based on the mean systemic arterial pressure; the second mechanism adapts the flow resistance in the pulmonary arterial system based on the mean pulmonary arterial pressure. Finally, the unstressed volume in the systemic veins is adapted based on the total cardiac output, i.e., the venous tone is adapted. The simulated inflow cannula of the tVAD was connected with the left ventricle, and the simulated outflow cannula was attached to the aorta. This model was used to generate the constant speed and the sinusoidal speed profile solutions by forward simulation, as described in Section 2.1.

Figure 1. The subsystems of the numerical circulation model described in [24,25], including a ventricular assist device. The optimal-control problem (OCP) considered in the current paper is based on the reduced model indicated by the gray box. There are three boundary conditions that were computed by the full model: the pulmonary arterial pressure, the systemic venous pressure and the systemic arterial resistance. The latter is influenced by the baroreflex mechanism. These quantities do not change, as long as the total cardiac output and, thus, the mean arterial pressure are unchanged.

In order to reduce the size and the complexity of the OCP described in Section 2.2, a reduced version of the model was used, as indicated by the gray box in Figure 1. Different choices for the boundaries of the reduced model are possible. Since the cardiac output was fixed and only the left ventricle was supported, the influence of the pump control strategy on the right ventricle was negligible. The boundaries need to be chosen, such that the pump speed profile does not have an effect on the pressure waveform at the boundaries. Choosing the pulmonary arterial pressure and the systemic venous pressure has been shown to fulfill these requirements. For these boundary conditions, the trajectories resulting from the constant speed solution were applied. Furthermore, the systemic arterial resistance from the same solution was enforced during the optimization. This setup ensured that the demanded cardiac output of 5Lmin 5Lmin led to the same mean systemic arterial pressure of 100mmHg 100mmHg, as in the solutions obtained from the full model.

2.1. Simulation of the Reference Solutions

The solutions for the constant speed and the sinusoidal speed profile cases were generated by forward simulations of the full model. Thereby, the mean pump speed, _{ωc ωc} , was regulated automatically, to yield the desired cardiac output of 5Lmin 5Lmin . In the constant speed case, there is one unique pump speed of _ωc=4,177rpm _ωc=4,177rpmthat provides the desired cardiac output.

The sinusoidal speed profiles are illustrated in Figure 2. They are defined as:

ω(t)=_ωc+_ωA·sin2·π·_γc(t)+φ+0.25

where _{ωA ωA} is the amplitude in rpm rpm and φ is the phase shift as a fraction of the duration of one cardiac cycle. The variable, _γc(t) _γc(t) , is a periodic, sawtooth-shaped signal representing the progress through each cardiac cycle. The variables, φ and _γc(t) _γc(t) , range from 0 to 1, where _γc(t)=0 _γc(t)=0 represents the onset of the ventricular systole. The phase offset of 0.25 0.25 was chosen, such that at a phase shift of φ=0 φ=0, the maximum tVAD speed occurs exactly at the onset of the ventricular systole.

Figure 2. Timing of the sinusoidal speed profile with respect to the cardiac cycle. The start of a cardiac cycle was defined as the onset of the ventricular systole (dotted vertical lines). At these instants, the cardiac cycle signal, γc(t) γc(t) , which was used to generate the sinusoidal speed profile, is reset. For the sinusoidal speed profiles, the two parameters, φ (phase shift) and ωA ωA (amplitude), are defined.

A “brute-force” analysis was performed: The two parameters, _{ωA ωA} and φ, were varied in small steps throughout their admissible range, and for each combination, a forward simulation was executed. The steps chosen were _ωA={250,500,…,3250,3500}rpm _ωA={250,500,…,3250,3500}rpm and φ={0,0.05,…,0.9,0.95} φ={0,0.05,…,0.9,0.95} . The data of the resulting 280 simulations were stored and can be used to evaluate any objective function, e.g., the one described by Equation (6a) below.

2.2. Optimal-Control Problem

This section is divided into four parts. Section 2.2.1 details the mathematical model of the reduced cardiovascular system. Section 2.2.2 defines the mathematical model of the blood pump, and Section 2.2.3 describes the formulation of the continuous-time OCP. Finally, Section 2.2.4 describes the transformation of this OCP into a nonlinear program (NLP), as well as the method for solving it.

2.2.1. Reduced Model of the Cardiovascular System

The differential equations describing the cardiovascular system are:

_x˙1(t)=_ppa(t)-_ppvx1(t)_Rpa-_ppvx1(t)-_plax2(t),t_Rpv

_x˙2(t)=_ppvx1(t)-_plax2(t),t_Rpv-_x3(t)

_x˙3(t)=1_Lmvplax2(t),t-_plvx(t),t+_ps,mv(t)-_Rmv·_x3(t)

_x˙4(t)=_x3(t)-_x5(t)-_x9(t)

_x˙5(t)=1_Lavplvx(t),t-_paox(t)+_ps,av(t)-_Rav·_x5(t)

_x˙6(t)=_x5(t)+_x9(t)-_x7(t)

_x˙7(t)=_paox(t)-_psax8(t)_Lsa

_x˙8(t)=_x7(t)-_psax8(t)-_psv(t)_Rsa

where each state variable describes either a volume or a volume flow in the cardiovascular system; see Table A1 in the Appendix. The right-hand sides of most of these differential equations depend on pressures in the cardiovascular system that are static functions of the state variables. The pressures in the pulmonary vein and in the systemic arteries are calculated assuming constant compliance and a volume offset:

_ppvx1(t)=_x1(t)-_V0,pvCpv

_psax8(t)=_x8(t)-_V0,saCsa

For the aortic pressure, visco-elasticity is taken into account, as well:

_paox(t)=_x6(t)-_V0,aoCao+_Rsa,i·_x5(t)-_x7(t)+_x9(t)

The pressures in the left atrium and the left ventricle are calculated using empirical time-varying elastance laws. For the left ventricle, a visco-elastic effect was included. Thus, the equations read:

_plax2(t),t=_Po,la+_φp,lax2(t),t+_ela(t)·_φa,lax2(t),t-_φp,lax2(t),t

_plvx(t),t=_Po,lv+_φp,lvx4(t),t+S·_elv(t)·_φa,lvx4(t),t-_φp,lvx4(t),t+_Ri,lv(t)·_x3(t)-_x5(t)-_x9(t)

The pressure-volume relations of the contracting chambers (left ventricle and left atrium) are defined as:

_φp,iVi(t)=_Emin,i·_Vi(t)-_V0,i+_KiVsat,i-_Vi(t)-_V0,i-_{Ki Vsat,i}

_φa,iVi(t)=1-^{_Vpmax,i-_Vi(t)_Vpmax,i-_V0,i2}·_Pmax,i

where the index, i, represents either the left atrium or the left ventricle. Accordingly, _Vi(t) _Vi(t) represents either the left atrial volume, _x2(t) _x2(t) , or the left ventricular volume, _x4(t) _x4(t) . The subscript, p p , in Equation (2n) denotes that the equation describes the passive pressure-volume relationship (relaxed muscle); the subscript, a a , in Equation (2o) denotes the active pressure-volume relationship (contracting muscle).

The model depends on several time-varying parameters: the pressures at the boundaries, _ppa(t) _ppa(t) and _psv(t) _psv(t) , the normalized elastance functions, _ela(t) _ela(t) and _elv(t) _elv(t) , and the internal resistance of the left ventricle, _{Ri,lv Ri,lv} . Figure 2 shows the two normalized elastance functions:

_ei(t)=12+12·cosπ-π_g1,i·_γc(t-_Ti),0≤_γc(t-_Ti)≤_g1,i12+12·cosπ·_g1,ig1,i-_g2,i-π_g1,i-_g2,i·_γc(t-_Ti), _g1,i<_γc(t-_Ti)≤_g2,i0,else

The index, i, denotes either the left ventricle (lv) or the left atrium (la). In the former case, the parameter, _{Tlv Tlv} , is equal to 0, and in the latter case the value of _Tla=0.2s _Tla=0.2s accounts for the time delay between the contraction of the two chambers. Table A2 in the Appendix shows the parameter values used in Equation (2).

The slack variables, _{ps,mv ps,mv} and _{ps,av ps,av} , were used to model the heart valves, which impose a strictly non-negative flow. To achieve this non-smooth behavior, the following set of constraints needs to be satisfied [28]:

_x3(t)·_ps,mv(t)=0

_x5(t)·_ps,av(t)=0

_x3(t)≥0

_ps,mv(t)≥0

_x5(t)≥0

_ps,av(t)≥0

2.2.2. Model of the Blood Pump

Two differential equations were used to describe the acceleration of the fluid in a mixed-flow blood pump (Deltastream DP2, Medos Medizintechnik AG, Stolberg, Germany) and the acceleration of the rotor shaft:

_x˙9(t)=1_LVAD·H_x9(t),_x10(t)-_paox(t),t-_plvx(t),t

_x˙10(t)=1_ΘVAD·-T_x9(t),_x10(t)+_kVAD·I(t)

Here, _x9(t) _x9(t) is the flow rate through the pump and _x10(t) _x10(t) is the rotational speed of the pump. The variable, I(t) I(t) , is the pump current. The “pressure head”, H, and the “hydraulic torque”, T, of the pump were stored in two-dimensional maps obtained from static measurements on the hybrid mock circulation. The parameter, _{kVAD kVAD} , relates the produced torque to the electric current and was supplied by the motor manufacturer. The fluid inertance, _{LVAD LVAD} , and the rotor inertia, _{ΘVAD ΘVAD} , were identified using dynamic measurements. The variables describing the pump model are summarized in Table A3 in the Appendix.

2.2.3. Formulation of the Optimal-Control Problem

The reduced model has _nx=10 _nx=10 state variables. The first eight ( _x1,...,_{x8 x1},...,_x8 ) describe the circulation Equation (2), whereas the two remaining ones ( _x9,_{x10 x9},_x10 ) describe the mixed-flow blood pump Equation (4). The model has one control input, namely the pump current, I(t) I(t) . The two slack variables, _{ps,mv ps,mv} and _{ps,av ps,av} , represent, from a structural point of view, two additional inputs. The control input and these two slack variables were gathered in the input vector, u(t) u(t) , which, consequently, has the dimension _nu=3 _nu=3 . All state variables were gathered in the state vector, x(t) x(t). The corresponding right-hand sides of the differential equations describing the model were stacked in the vector-valued “model function”, f:

x˙(t)=fx(t),u(t),t

The mathematical representation of the optimization problem to be solved comprises an objective function Equation (6a) and several constraints Equations (6b)–(6i):

min.x(t),u(t)Jx(t),t=_∫tbtf Lx(t),tdt

s.t.x˙(t)=fx(t),u(t),t,

sx(t)=0

gx(t)=0

x(_tb)=x(_tf)

u(_tb)=u(_tf)

_xlb≤x(t)≤_xub

_ulb≤u(t)≤_uub

_tb≤t≤_tf.

The function, L, in Equation (6a) is an arbitrary function of the state variables and is described below. The time instances, _{tb tb} and _{tf tf} , denote the start of two consecutive cardiac cycles, where the heart rate is equal to 90bpm 90bpm , i.e., _tf-_tb=0.667s _tf-_tb=0.667s . The vector-valued continuity constraint Equation (6b) ensures that, for a solution of the OCP, the model equations are satisfied. The constraint Equation (6c) represents Equations (3a) and (3b), which model the heart valves. The constraint Equation (6d) ensures the desired cardiac output of _V*CO=5Lmin _V*CO=5Lmin,

gx(t)=1_tf-_{tb∫tbtf x5}(t)+_x9(t)dt-_V*CO

A periodic solution was sought, defined by the constraints Equations (6e) and (6f). The constraints Equations (6g) and (6h) implement the upper and lower bounds on the state variables and the inputs, as specified in Table A1 in the Appendix. Different values were chosen for the lower bound on the pump flow, _{x9,lb x9,lb} . The standard value of 0 implies that no backflow occurs from the aorta to the left ventricle. The influence on the solution when a strictly positive flow was enforced is analyzed in Section 3.4.

The objective function Equation (6a) was chosen in accordance with physiological considerations. The application of a VAD can lead to a permanent closure of the aortic valve. Fusion of the aortic valvular cusps and a thrombus formation at the aortic valve may occur in this case [29,30,31]. Including the flow through the aortic valve as the first component of the objective function allows solutions to be enforced where the flow through the aortic valve is maximized and, thus, opening of the valve is ensured. A further goal of the installation of a VAD is the unloading of the heart. Thus, a reduction of the hydraulic work the left ventricle needs to deliver was chosen as a second component of the objective function. The main task of a VAD is to maintain the perfusion of the patient, which was enforced by the perfusion constraint Equation (6d).

The tradeoff between the two conflicting components of the objective function is adjusted by the weighting parameter, ρ∈[0,1] ρ∈[0,1]. The objective function thus reads:

_L1x5(t)=-_x5(t)·_L1,N

_L2x(t),t=_plvx(t),t·_x5(t)+_x9(t)-_x3(t)·_L2,N

Jx(t),t,ρ=_∫tbtf1-ρ·_L1x5(t)+ρ·_L2x(t),tdt

The parameter _L1,N=^10-2J·smL _L1,N=^10-2J·smL normalizes _{L1 L1} , such that the typical values of both components of the objective function are within the same numerical range. The parameter _{L2,N L2,N} =133.3 × ^{10-6 10-6} converts the units of _{L2 L2} from mmHg·mL mmHg·mLto J.

2.2.4. Numerical Solution of the Optimal-Control Problem

Direct transcription was applied to solve the OCP Equation (6), transforming the continuous-time problem into an NLP [32,33]. Despite its large size, the resulting NLP can be solved efficiently, due to its structure and the performance of current NLP solvers. Furthermore, this approach provides a straightforward handling of any type of constraints [34].

The family of Radau collocation schemes was applied in this work [35]. Collocation methods represent each state variable _xi(t) _xi(t) as a polynomial and, thus, provide a continuous solution. When used in the context of direct transcription, their implicit nature is irrelevant. “Radau IIA” methods were chosen because to their advantageous stability properties [36]. For the problem at hand, using a low-order method with a fine time discretization yielded the best results concerning the convergence of the NLP solver and the avoidance of local minima. In fact, the results shown below were obtained using the first-order representative of this collocation family, which is the well-known Euler-backward scheme.

All of the equations forming the OCP Equation (6) have to be discretized consistently. A grid of N points was used, yielding N-1 N-1 intervals of length _hk=_tk+1-_{tk hk}=_tk+1-_tk . The notation _xk=x(_tk) _xk=x(_tk) was adopted. The continuity constraint Equation (6b) was transformed to discrete defect constraints; all integrals became weighted sums, and the constraints on the inputs and on the state variables were imposed at the discretization points. Therefore, the continuous OCP Equation (6) is transformed to:

min._x•,_u•Jd=∑k=1N-1_hk-1·L(_xk,_tk)

whereL(_xk,_tk)=(ρ-1)·_x5,k·_L1,N+ρ·_plv(_xk,_tk)·(_x5,k+_x9,k-_x3,k)·_L2,N

s.t._xk+1=_xk+_hk·f(_xk+1,_uk+1),k=0…N-2

_sk=_x3,k·_u2,kx5,k·_u3,k=0,k=0…N-1

∑k=1N-1_hk-1·(_x5,k+_x9,k)=(_tN-1-_t0)·_V*CO

_x0=_xN-1

_u0=_uN-1

_xlb≤_xk≤_xub,k=0…N-1

_ulb≤_uk≤_uub,k=0…N-1

_t0=_tb,_t1,…,_tN-2,_tN-1=_tf

Equation (9) defines an NLP in n=N·(_nu+_nx) n=N·(_nu+_nx) variables, namely _{uk uk} and _{xk xk} for k=0…N-1 k=0…N-1 , which has _nc=(N-1)·_nx+2·N+1+(_nx+_nu) _nc=(N-1)·_nx+2·N+1+(_nx+_nu) equality constraints. Its objective function Equation (9a) is nonlinear in the NLP variables, and it has nonlinear constraints Equations (9b) and (9c), linear constraints Equations (9d), (9e) and (9f), as well as simple bounds Equations (9g) and (9h). All state variables and inputs were scaled, such that the range between the lower and the upper bounds was mapped to the interval, [0,1] [0,1]. This normalization ensures a similar magnitude of all NLP variables. The resulting improvement in numerical condition of the problem increased the reliability and the convergence speed of the NLP solver.

The equality-constrained NLP is commonly solved using a Newton-type iteration. Due to the required factorizations (or re-factorizations, respectively) of the large matrices involved, a typical runtime between O(ⁿ²) O(ⁿ²) and O(ⁿ³) O(ⁿ³) results [37]. The second-order derivatives were approximated by an iterative update using first-order derivative information. Two factors are thus key for a fast solution of the OCP: (a) the reduction of the problem size and (b) the exploitation of the problem structure to reduce the number of evaluations of the objective and model functions. The following paragraphs outline the procedure to achieve both.

The problem size was kept small by first using a discretization with a constant step size of h=20ms h=20ms , resulting in N=34 N=34 discretization points. The resulting solution was then used to initialize a refined version of the problem. Namely, the step size was halved, which results in an NLP that is twice as large. However, since a close initial guess is available, only a few iterations are required to solve this NLP to the demanded tolerance. For all solutions presented below, two refinements were performed, which yielded a final step size of 5ms 5ms.

The Jacobian matrix contains the partial derivatives of the objective function Equation (9a) and of all constraints Equations (9b)–(9f) with respect to the NLP variables. Since the objective function and each constraint depend on a few NLP variables only, most of its entries are zero. This fact is termed “sparsity”, and its exploitation is the key to an efficient numerical solution of directly transcribed OCPs [38]. Most importantly, the objective function, L, and the model function, f, only occur in Equations (9a) and (9b), respectively. Furthermore, their arguments are the state variables and the inputs at the same discretization point. Therefore, all derivative information required for the construction of the Jacobian matrix is provided by the partial derivatives of the model function and of the objective function with respect to the state variables, as well as the inputs at each discretization node.

Any NLP solver can calculate the Jacobian matrix by numerical differentiation. However, depending on its algorithmic implementation, it is unable to fully exploit the problem sparsity. In order to always obtain the least possible number of evaluations of the objective and the model functions, the Jacobian matrix was constructed by custom code and was then passed to the NLP solver. Linear terms Equations (9d)–(9f) have constant derivatives, whereas simple bounds Equations (9g)–(9h) are handled directly by the NLP solver. The derivatives of objective function Equation (9a) and of defect constraints Equation (9b) were constructed from the derivatives of the objective and the model functions, which were calculated analytically.

As an NLP solver, SNOPT was used [39]. Iteratively refining the problem discretization and providing the Jacobian matrix to the solver reduced the typical runtime for solving one OCP from about 30min 30min to less than 1min 1min. This reduction in computational time enabled large-scale parametric studies and a quick comparison of various problem formulations.

2.3. Validation on the Hybrid Mock Circulation

The hybrid mock circulation presented in [25] was used to validate one solution resulting from the optimization procedure described in Section 2.2.4. This test bench interfaces a hydraulic system with the full model that is used for the simulations. The model runs in real time on a PC. The calculated left ventricular and aortic pressures are reproduced by two pressure-controlled hydraulic reservoirs. The mixed-flow blood pump is connected to these reservoirs and, thus, perceives pressures that closely approximate an in vivo environment. The flow rate through the blood pump is measured using an ultrasonic flow probe (11PXL, Transonic Systems Inc., Ithaca, NY, USA). It is fed back to the circulation model to allow a real-time interaction between the blood pump and the numerical circulation model.

An optimized speed trajectory was realized using a modified speed controller for the blood pump. After several cardiac cycles, stationary operation was reached, as imposed within the OCP by periodicity conditions Equations (6e) and (6f). Upon reaching steady state, the data from the last cardiac cycle was stored for further analysis. The pressures in the two reservoirs, representing the up- and down-stream pressures of the blood pump, as well as the pump speed and flow rate were measured. In addition, all relevant data from the numerical simulation were saved along the measured signals.

3. Results

This section is organized as follows. Section 3.1 describes the runtime and the convergence of the numerical solution of the OCP Equation (9). Section 3.2 presents the measurement results of one solution of the OCP applied on the hybrid mock circulation in order to validate the model of the blood pump and the simplified model of the blood circulation. Section 3.3 presents the results of the sinusoidal speed profile investigation. Finally, in Section 3.4, both the constant speed solution and the optimal sinusoidal solution are compared to speed profiles obtained by the numerical solution of OCP Equation (9) for several choices of the tradeoff parameter, ρ, and the lower bound on the pump flow, _{x9,lb x9,lb}.

3.1. Runtime and Convergence of the Optimization

Figure 3 illustrates the convergence properties of the method used. Table 1 lists the content of panels (a) and (c) of Figure 3 in numerical form. The NLP Equation (9) was solved with ρ={0,0.5,1} ρ={0,0.5,1} . Solutions were calculated for different uniform step sizes h≈{85,30,15,7,4,3}ms h≈{85,30,15,7,4,3}ms . The corresponding numbers of NLP variables are n={150,350,600,1500,2300,3000} n={150,350,600,1500,2300,3000} . Figure 3(a) shows the root mean squared difference between two successive solutions of the state variable, _x9(t) _x9(t) , which is the flow rate through the blood pump. Figure 3(c) shows the corresponding runtimes the NLP solver needed to converge. Panels (c) and (d) show three solutions with 150, 600 and 3,000 3,000 NLP variables, respectively, for ρ=1 ρ=1 and ρ=0 ρ=0 . Clearly, the obtained solution did only change marginally when a problem size of n≈1,000 n≈1,000was exceeded. A solution to a problem of this size can be obtained within less than one minute.

3.2. Validation of the Model Reduction

Figure 4 shows the graphs resulting from the solution of the OCP with ρ=0 ρ=0 (in silico) and the reproduction of this optimal solution on the hybrid mock circulation (in vitro). The current, the speed and the flow rate of the pump, as well as the relevant pressures are shown. In addition, the flow through the aortic valve and the left ventricular pressure-volume diagrams are displayed. When the optimal solution was reproduced on the hybrid mock circulation, the total cardiac output amounted to 4.92Lmin 4.92Lmin . Compared to the demanded value of 5Lmin 5Lmin , which was exactly matched by the in silico solution, this deviation of -1.6% -1.6%is negligible.

In the optimal solution, the pump actuation was such that the pump flow rate was close to zero during ventricular systole. Therefore, the pressure-volume diagram contains a nearly isovolumetric contraction phase followed by the ejection of blood through the aortic valve. The current, the speed and the flow rate of the pump all reached their lower or upper bounds at certain instants in time; see Figure 4.

The boundary conditions used for solving the OCP are the pulmonary arterial pressure, the systemic venous pressure and the systemic arterial resistance, as indicated in Figure 1. The trajectories from the simulation using a constant pump speed were compared with the signals measured when the optimal solution was reproduced on the hybrid mock circulation. The maximum absolute difference in the pulmonary arterial pressure amounted to 0.29mmHg 0.29mmHg and the root mean squared error was equal to 0.04mmHg 0.04mmHg . The respective values for the systemic venous pressure were 0.04mmHg 0.04mmHg and 9.5 × ^{10-4 10-4} mmHg mmHg . The systemic arterial resistance imposed by the baroreflex mechanism is assumed to be constant for the OCP. Its value changed from 1.113mmHg·smL 1.113mmHg·smL for the constant pump-speed simulation to 1.110mmHg·smL 1.110mmHg·smL for the in vitro reproduction of the optimal solution. This difference of just -0.27% -0.27% indicates that the use of the reduced model within the OCP is valid. The unavoidable modeling errors in Equation (4), as well as the tracking error of the pressure control system at the hybrid mock circulation led to differences between the optimized solution and the measured data, as shown in Figure 4. These effects accumulated to a reduction of the end-diastolic left ventricular volume from 133mL 133mL to 128mL 128mL and a reduction of the flow through the aortic valve from 1.08Lmin 1.08Lmin to 0.88Lmin 0.88Lmin.

Figure 3. Convergence of the nonlinear program (NLP) solver for NLP Equation (9) and the resulting pump flow profiles. Panel (a) shows the RMS of the difference between a solution and the next finer discretized solution. Panel (b) shows the pump flow trajectories obtained for ρ=1 ρ=1 and three differently accurate transcriptions. Panel (c) shows the runtime the NLP solver needs to converge, depending on the problem size. Panel (d) shows the same information as panel (b) for ρ=0 ρ=0 .

[ Table omitted. See PDF. ]

Figure 4. The optimal solution for maximizing the flow through the aortic valve ( ρ=0 ρ=0 ), in silico and in vitro. Panels (a)–(c) show the current, the speed, and the flow rate of the pump. Panel (d) shows the left ventricular (solid) and the aortic (dashed) pressures. Panel (e) shows the flow through the aortic valve. Panel (f) shows the left ventricular pressure-volume diagrams.

3.3. Analysis of the Sinusoidal Speed Profiles

Figure 5 shows the results from the parametric study described in Section 2.1. Contour lines of the mean pump speed, the minimum pump flow rate, the mean flow through the aortic valve and the left ventricular stroke work are plotted over the parameter space _ωA×φ _ωA×φ . As stated above, the mean pump speed was adjusted, such that the demanded cardiac output of 5Lmin 5Lminis provided.

Figure 5. Results of the parametric study on the sinusoidal speed profile. All plots show contour lines of a specific measure, as indicated in the corresponding heading. Panel (a) shows the mean pump speed, ωc ωc in rpm rpm . Panel (b) shows the minimum value of the flow rate through the pump in Lmin Lmin . The gray area indicates the admissible region. Backflow from the aorta to the left ventricle occurs in the complementary region. Panel (c) shows the mean flow rate through the aortic valve in Lmin Lmin . There is a large region where the aortic valve remains closed. Panel (d) shows the left ventricular stroke work. The asterisks and arrows in panels (c) and (d) denote the optimal point according to objective function Equation (9a) and the corresponding directions of the steepest descent. The admissible region shown in panel (b) is repeated in panels (c) and (d).

The mean pump speed varied between 3,400rpm 3,400rpm and 4,700rpm 4,700rpm in order to provide the desired cardiac output. The minimum pump flow rate decreased with increasing amplitude, which caused an increased backflow from the aorta to the left ventricle. Combined with the requirement that there must be no backflow through the pump, i.e., _x9,lb=0 _x9,lb=0 , this relationship restricted the admissible range of the amplitude to the gray area shown in Figure 5(b). For positive choices of the parameter, _{x9,lb x9,lb} , this area became even smaller. Clearly, there is a large region in the parameter space where the aortic valve remains closed, as shown in Figure 5(c). Finally, the points of minimum left ventricular stroke work and maximum flow through the aortic valve do not coincide.

For all pairs of amplitude and phase shift, the value of the objective function Equation (6a) was calculated from the corresponding simulation data. From all pairs that fulfill constraints Equations (6e)–(6h), the one with the minimum value of the objective function was selected as the optimal choice. This solution, defined by φ=0.8 φ=0.8 and _ωA=2,000rpm _ωA=2,000rpm , is indicated in Figure 5(c,d) by an asterisk. The structure of the system is such that the location of the optimal sinusoidal speed profile is invariant with respect to the tradeoff parameter, ρ. This finding stems from the fact that for both the flow through the aortic valve and the left ventricular stroke work, the direction of decreasing objective-function values (arrows) is towards the non-admissible region in which backflow through the pump occurs.

3.4. Comparison of the Numerically Optimized Speed Profiles with Reference Speed Profiles

The constant speed solution together with the unique optimal solution for the sinusoidal speed profile is presented in Figure 6. Alongside, the results from the optimization procedure for ρ={0,0.5,1} ρ={0,0.5,1} are displayed. In contrast to the unique optimal sinusoidal speed profile, the individual solutions from the optimization differed for variations of the parameter, ρ. As a consequence of the shift of the tradeoff in the objective function, the stroke volume and the stroke work of the left ventricle varied, as shown in Figure 6(h,j).

Figure 7 illustrates the tradeoff between the mean flow through the aortic valve and the left ventricular stroke work. The values for the constant speed solution, the optimal sinusoidal speed profile and the OCP solutions are shown. For the latter two cases, three different choices for _{x9,lb x9,lb} were used. Clearly, compared to the standard case ( _x9,lb=0 _x9,lb=0 ), the values of the attainable objective-function got worse when stasis was avoided in the pump cannulas ( _x9,lb=1Lmin _x9,lb=1Lmin or _x9,lb=2Lmin _x9,lb=2Lmin ). The restriction of the admissible minimum pump flow indirectly also increases the minimum pump speed that resulted. For the case _x9,lb=1Lmin _x9,lb=1Lmin , the minimum pump speed varied between 0rpm 0rpm and 550rpm 550rpm , and for the case _x9,lb=2Lmin _x9,lb=2Lmin , the minimum pump speed varied between 570rpm 570rpm and 990rpm 990rpm. In each case, the optimal solutions performed better than the sinusoidal speed profiles in terms of minimizing the value of the objective function. The performance of the constant-speed solution clearly was the least satisfactory.

4. Discussion

In the current study, numerical optimal control was used to optimize the functional interaction between tVADs and the circulatory system. The performance of speed profiles was assessed by a physiologically motivated objective function. A mathematical model was used during the optimization, and the results were validated on a real blood pump in a hybrid mock circulation. This approach yields a better performance than constant-speed or sinusoidal profiles.

The necessity of taking into account multiple physiological constraints simultaneously has been pointed out in previous work [14]. In that study, bounds were specified for the stroke volume of the left ventricle, for the mean left atrial and ventricular pressures, as well as for the arterial pulse pressure. All results from the optimal-control approach presented in the current paper stay within these bounds. The framework of numerical optimal control can handle any type of constraints imposed on the state variables and the inputs. Whenever certain physiologic measures appear to be out of range, additional constraints can be specified in the OCP. Such constraints would either force the measures in question to stay within the desired limits, or, if no feasible solution is found, an improper selection of these limits would be revealed. This feature highlights the benefits of the methodology presented, inasmuch as it represents a flexible framework where the objective function and the constraints can be adapted to specific requirements.

Figure 6. Panels (a–e) show the pump speed, the pump flow rate, the flow through the aortic valve, the left ventricular pressure and the left ventricular pressure-volume diagram for the constant speed solution and for the optimal sinusoidal speed profile. The latter is characterized by φ=0.8 φ=0.8 , ωA=2,000rpm ωA=2,000rpm and ωc=3,968rpm ωc=3,968rpm . Panels (f–j) show the trajectories of these quantities for the solutions of the OCP with ρ={0,0.5,1} ρ={0,0.5,1} . For all cases, a value of x9,lb=0Lmin x9,lb=0Lmin is used. The optimized solution with ρ=0 ρ=0 is equal to the optimized solution (in silico) shown in Figure 4.

Figure 7. Tradeoff between the flow through the aortic valve and the left-ventricular stroke work. Recall that a higher flow through the aortic valve and a lower stroke work decrease the objective-function value Equation (6a). Values for the constant speed solution, the optimal sinusoidal speed profiles and the optimal solutions are shown. Three different values for the lower bound on the minimum pump flow rate are applied. The parameter, ρ, was varied between zero (the flow through the aortic valve is maximized) and one (the left ventricular stroke work is minimized). The optimal parameter set for the sinusoidal profile is invariant w.r.t.ρ; see Figure 5.

As described in the introduction, numerous approaches exist for speed modulation of tVADs, synchronized or asynchronous to the ventricular contraction and with different types of speed profiles. None of these approaches has been applied clinically as yet, where, normally, a constant speed is applied. The bandwidth of the reference-speed tracking in a real application is important. For example, square-wave speed profiles cannot be tracked exactly, due to the rotor inertia and the limited electric power that is available. When speed profiles are derived from the solution of an OCP, these limitations can be taken into account by the introduction of appropriate constraints.

The OCP Equation (6) is formulated for one cardiac cycle, i.e., it is assumed that the optimal solutions are periodic with the same duration as the cardiac cycle. To verify this assumption, an OCP similar to Equation (6) was solved, where the time horizon included two cardiac cycles. No periodicity constraints were imposed in between these two cycles. The resulting solution was periodic with the duration of only one cardiac cycle, and the speed profiles for the two cardiac cycles were identical.

When such synchronized speed profiles are to be applied in clinical practice, the timing of the ventricular contraction must be extracted from a measured signal, which might require data from an additional sensor, such as from the electrocardiogram [40]. Similar methods were available in some first-generation, volume displacement VADs [41]. Potentially, the synchronization of a VAD to the cardiac cycle could be realized using pump-intrinsic signals only by extending the work presented in [42]. Nevertheless, when speed modulation is to be applied clinically, the option of using optimized speed profiles as presented in the current paper should be considered, due to its superior performance, as indicated in Figure 7.

The solutions generated by the optimal-control approach are non-causal and rely on an offline, model-based optimization of the pump-speed profile. Two options for realizing their potential in the context of causal control, i.e., in a clinical application with the possibility for patient-specific treatment, are outlined next.

(1) The optimal speed profiles can be parametrized for various physiological states, e.g., various heart rates or changing objective functions. These pre-calculated optimal solutions could then be programmed in a clinical controller. A possible application would be to switch between different speed profiles, including the standard solution where a constant speed is applied. Further research is necessary to verify this approach.

(2) Online algorithms could be developed on the basis of the insights gained from the optimal solutions. For example, Figure 6 reveals that for minimizing the objective function Equation (6a), a two-step control strategy should be applied: In the first part of the speed profile, the pump speed is chosen, such that the pump flow stays at its lower boundary. In the second part of the speed profile, the pump flow needs to be controlled, such that the desired cardiac output is reached and, therefore, a flow pulse is generated. The parameter, ρ, thereby adopts the role of a phase shift and varies the timing of the two parts with respect to the ventricular contraction.

Furthermore, the solutions generated by solving an OCP can be used for benchmarking other approaches based on mathematical models. The performance of any control strategy can be compared to the optimal value, which allows one to assess the remaining potential. Finally, a fair comparison among various VAD types can be performed by comparing the optimal behavior of the devices in question.

The optimal pump flow profiles obtained in the current paper are comparable to those typically obtained with a pulsatile volume-displacement VAD, although with a lower amplitude. The effect of the timing of the ejection phase of such a pulsatile volume-displacement VAD on the stroke work and the aortic valve flow was demonstrated earlier in vivo [43]. Similarly to the results in the current paper, that analysis shows that a change in phase shift is required when the maximization of the aortic valve flow was assigned preference over the minimization of the left ventricular stroke work.

The blood damage induced by blood pumps represents a complex topic, and different outcomes in this regard have been observed clinically with different types of devices [44]. Increased hemolysis due to speed modulation has been reported for a commercially available centrifugal pump in vitro [45,46]. Another group designed an impeller to reduce blood damage during speed modulation, which was successfully tested in vivo [47,48]. A recently designed VAD (HeartMate III, Thoratec Corp., Pleasanton, CA, USA) is expected to go into clinical trial soon and employs speed modulation with gradients of up to 24,000rpms 24,000rpms [49]. So far, no mathematical description for blood damage by means of lumped parameter models exists. Therefore, such a description could not be included into the framework presented as yet.

The results shown in the current manuscript demonstrate the potential of the methodology proposed. Future studies should focus on a sensitivity analysis of the optimal solutions. On the one hand, the benefit attainable by the application of optimized speed profiles has to be assessed on a more comprehensive range of the state of the cardiovascular system. On the other hand, the sensitivity of the objective-function value to deviations from a design point used during the optimization should be analyzed. This analysis indicates which physiological parameters need to be included as scheduling parameters in a clinical application. At the same time, the trends of the optimal operating strategy w.r.t. these parameters are revealed, from which guidelines for the development of online control strategies may be derived.

The method presented is able to generate optimal speed profiles for tVADs. By integrating the cardiovascular system into the optimization procedure, this approach not only allows the operation of the tVAD itself to be optimized, but it also improves its interaction with the cardiovascular system. The methodology presented is well suited as a tool for further research of speed modulation for tVADs. Transferring the results obtained to the clinical setting would allow one to adjust the pump performance according to the intention-to-treat of the individual patient, such as the recovery of the native heart or lifelong support.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix

Table A1 provides an overview of the state variables, the control inputs and the slack variables used in the reduced circulation model and the pump model. Table A2 lists the parameters used in reduced numerical circulation model Equation (2). Table A3 details the parameters of the model of blood pump Equation (4).

Table A1. Description of the state variables and the inputs of the reduced circulation model and the pump, including their lower (LB) and upper (UB) bounds.

[ Table omitted. See PDF. ]

^*_x9,lb={0,1,2}Lmin={0,16.6,33.3}mLs ^*_x9,lb={0,1,2}Lmin={0,16.6,33.3}mLs

Table A2. Parameters of the reduced circulation model.

[ Table omitted. See PDF. ]

Table A3. Parameters of the pump model.

[ Table omitted. See PDF. ]

References

1. Timms, D. A review of clinical ventricular assist devices. Med. Eng. Phys. 2011, 33, 1041–1047.

2. Burkhoff, D.; Klotz, S.; Mancini, D.M. LVAD-induced reverse remodeling: Basic and clinical implications for myocardial recovery. J. Card. Fail. 2006, 12, 227–239.

3. Drakos, S.G.; Terrovitis, J.V.; Anastasiou-Nana, M.I.; Nanas, J.N. Reverse remodeling during long-term mechanical unloading of the left ventricle. J. Mol. Cell. Cardiol. 2007, 43, 231–242.

4. Kirklin, J.K.; Naftel, D.C.; Kormos, R.L.; Stevenson, L.W.; Pagani, F.D.; Miller, M.A.; Baldwin, J.T.; Young, J.B. Fifth INTERMACS annual report: Risk factor analysis from more than 6,000 mechanical circulatory support patients. J. Heart Lung Transpl. 2013, 32, 141–156.

5. Yoshizawa, M.; Takeda, H.; Watanabe, T.; Miura, M.; Yambe, T.; Katahira, Y.; Nitta, S. An automatic control algorithm for the optimal driving of the ventricular-assist device. IEEE Trans. Biomed. Eng. 1992, 39, 243–252.

6. Arai, H.; Fujiyoshi, K.; Sakamoto, T.; Suzuki, A.; Swartz, M.T. Optimal control algorithm for pneumatic ventricular assist devices: Its application to automatic control and monitoring of ventricular assist devices. Artif. Organs 1996, 20, 1034–1041.

7. Yu, Y.C.; Boston, J.R.; Simaan, M.A.; Antaki, J.F. Minimally invasive estimation of systemic vascular parameters for artificial heart control. Control Eng. Pract. 2002, 10, 277–285.

8. Schima, H.; Vollkron, M.; Jantsch, U.; Crevenna, R.; Roethy, W.; Benkowski, R.; Morello, G.; Quittan, M.; Hiesmayr, M.; Wieselthaler, G. First clinical experience with an automatic control system for rotary blood pumps during ergometry and right-heart catheterization. J. Heart Lung Transpl. 2006, 25, 167–173.

9. Wu, Y.; Allaire, P.E.; Tao, G.; Olsen, D. Modeling, estimation, and control of human circulatory system with a left ventricular assist device. IEEE Trans. Control Syst. Technol. 2007, 15, 754–767.

10. Ferreira, A.; Boston, J.R.; Antaki, J.F. A control system for rotary blood pumps based on suction detection. IEEE Trans. Biomed. Eng. 2009, 56, 656–665.

11. Bearnson, G.B.; Olsen, D.B.; Khanwilkar, P.S.; Long, J.W.; Allaire, P.E.; Maslen, E.H. Pulsatile operation of a centrifugal ventricular assist device with magnetic bearings. ASAIO J. 1996, 42, M620–M624.

12. Bourque, K.; Dague, C.; Farrar, D.; Harms, K.; Tamez, D.; Cohn, W.; Tuzun, E.; Poirier, V.; Frazier, O.H. In vivo assessment of a rotary left ventricular assist device-induced artificial pulse in the proximal and distal aorta. Artif. Organs 2006, 30, 638–642.

13. Umeki, A.; Nishimura, T.; Takewa, Y.; Ando, M.; Arakawa, M.; Kishimoto, Y.; Tsukiya, T.; Mizuno, T.; Kyo, S.; Ono, M.; et al. Change in myocardial oxygen consumption employing continuous-flow LVAD with cardiac beat synchronizing system, in acute ischemic heart failure models. J. Artif. Organs 2013, 16, 119–128.

14. Shi, Y.; Brown, A.G.; Lawford, P.V.; Arndt, A.; Nuesser, P.; Hose, D.R. Computational modelling and evaluation of cardiovascular response under pulsatile impeller pump support. Interface Focus 2011, 1, 320–337.

15. He, P.; Bai, J.; Xia, D.D. Optimum control of the Hemopump as a left-ventricular assist device. Med. Biol. Eng. Comput. 2005, 43, 136–141.

16. Cox, L.G.E.; Loerakker, S.; Rutten, M.C.M.; de Mol, B.A.J.M.; van de Vosse, F.N. A mathematical model to evaluate control strategies for mechanical circulatory support. Artif. Organs 2009, 33, 593–603.

17. Shiose, A.; Nowak, K.; Horvath, D.J.; Massiello, A.L.; Golding, L.A.R.; Fukamachi, K. Speed modulation of the continuous-flow total artificial heart to simulate a physiologic arterial pressure waveform. ASAIO J. 2010, 56, 403–409.

18. Vandenberghe, S.; Segers, P.; Antaki, J.F.; Meyns, B.; Verdonck, P.R. Hemodynamic modes of ventricular assist with a rotary blood pump: Continuous, pulsatile, and failure. ASAIO J. 2005, 51, 711–718.

19. Pirbodaghi, T.; Weber, A.; Axiak, S.; Carrel, T.; Vandenberghe, S. Asymmetric speed modulation of a rotary blood pump affects ventricular unloading. Eur. J. Cardio-Thorac. 2013, 43, 383–388.

20. Pirbodaghi, T.; Axiak, S.; Weber, A.; Gempp, T.; Vandenberghe, S. Pulsatile control of rotary blood pumps: Does the modulation waveform matter? J. Thorac. Cardiovasc. Surg. 2012, 144, 970–977.

21. Tavoularis, S.; Sahrapour, A.; Ahmed, N.U.; Madrane, A.; Vaillancourt, R. Towards optimal control of blood flow in artificial hearts. Cardiovasc. Eng. 2003, 8, 20–31.

22. Tasch, U.; Geselowitz, D.B.; Sinha, A.; Hsu, H.K. A novel output feedback pusher plate controller for the Penn State electric ventricular assist device. J. Dyn. Syst. Meas. Control 1989, 111, 69–74.

23. Klute, G.K.; Tasch, U.; Geselowitz, D.B. An optimal controller for an electric ventricular-assist device: Theory, implementation, and testing. IEEE Trans. Biomed. Eng. 1992, 39, 394–403.

24. Colacino, F.M.; Moscato, F.; Piedimonte, F.; Arabia, M.; Danieli, G.A. Left ventricle load impedance control by apical VAD can help heart recovery and patient perfusion: A numerical study. ASAIO J. 2007, 53, 263–277.

25. Ochsner, G.; Amacher, R.; Amstutz, A.; Plass, A.; Schmid Daners, M.; Tevaearai, H.; Vandenberghe, S.; Wilhelm, M.J.; Guzzella, L. A novel interface for hybrid mock circulations to evaluate ventricular assist devices. IEEE Trans. Biomed. Eng. 2013, 60, 507–516.

26. Guyton, A.C.; Hall, J.E. Textbook of Medical Physiology, 12th ed.; Saunders: Philadelphia, PA, USA, 2010.

27. Suga, H.; Sagawa, K. Instantaneous pressure-volume relationships and their ratio in the excised, supported canine left ventricle. Circ. Res. 1974, 35, 117–126.

28. Pfeiffer, F. Deregularization of a smooth system—Example hydraulics. Nonlinear Dyn. 2007, 47, 219–233.

29. Rose, A.G.; Park, S.J.; Bank, A.J.; Miller, L.W. Partial aortic valve fusion induced by left ventricular assist device. Ann. Thorac. Surg. 2000, 70, 1270–1274.

30. John, R.; Mantz, K.; Eckman, P.; Rose, A.; May-Newman, K. Aortic valve pathophysiology during left ventricular assist device support. J. Heart Lung Transpl. 2010, 29, 1321–1329.

31. May-Newman, K.; Enriquez-Almaguer, L.; Posuwattanakul, P.; Dembitsky, W. Biomechanics of the aortic valve in the continuous flow VAD-assisted heart. ASAIO J. 2010, 56, 301–308.

32. Betts, J.T. Survey of numerical methods for trajectory optimization. J. Guid. Control Dynam. 1998, 21, 193–207.

33. Rao, A.V. A Survey of Numerical Methods for Optimal Control. In Proceedings of the AAS/AIAA Astrodynamics Specialist Conference, Pittsburgh, PA, USA, 9–13 August 2009.

34. Kameswaran, S.; Biegler, L.T. Simultaneous dynamic optimization strategies: Recent advances and challenges. Comput. Chem. Eng. 2006, 30, 1560–1575.

35. Hairer, E.; Wanner, G. Stiff differential equations solved by Radau methods. J. Comput. Appl. Math. 1999, 111, 93–111.

36. Butcher, J. Numerical Methods for Ordinary Differential Equations; John Wiley & Sons: Winchester, UK, 2003.

37. Kirches, C.; Bock, H.G.; Schloder, J.P.; Sager, S. A factorization with update procedures for a KKT matrix arising in direct optimal control. Math. Prog. Comput. 2011, 3, 319–348.

38. Patterson, M.A.; Rao, A.V. Exploiting sparsity in direct collocation pseudospectral methods for solving optimal control problems. J. Spacecr. Rockets 2012, 49, 364–377.

39. Gill, P.E.; Murray, W.; Saunders, M.A. SNOPT: An SQP algorithm for large-scale constrained optimization. SIAM Rev. 2005, 47, 99–131.

40. Amacher, R.; Ochsner, G.; Ferreira, A.; Vandenberghe, S.; Schmid Daners, M. A robust reference signal generator for synchronized ventricular assist devices. IEEE Trans. Biomed. Eng. 2013, 60, 2174–2183.

41. Farrar, D.J.; Compton, P.G.; Lawson, J.H.; Hershon, J.J.; Hill, J.D. Control modes of a clinical ventricular assist device. IEEE Eng. Med. Biol. 1986, 5, 19–25.

42. Moscato, F.; Granegger, M.; Edelmayer, M.; Zimpfer, D.; Schima, H. Continuous monitoring of cardiac rhythms in left ventricular assist device patients. Artif. Organs 2013, in press.

43. Amacher, R.; Weber, A.; Brinks, H.; Axiak, S.; Ferreira, A.; Guzzella, L.; Carrel, T.; Antaki, J.F.; Vandenberghe, S. Control of ventricular unloading using an electrocardiogram-synchronized Thoratec paracorporeal ventricular assist device. J. Thorac. Cardiovasc. Surg. 2013, 146, 710–717.

44. Heilmann, C.; Geisen, U.; Benk, C.; Berchtold-Herz, M.; Trummer, G.; Schlensak, C.; Zieger, B.; Beyersdorf, F. Haemolysis in patients with ventricular assist devices: Major differences between systems. Eur. J. Cardio-Thorac. 2009, 36, 580–584.

45. Tayama, E.; Niimi, Y.; Takami, Y.; Ohashi, Y.; Ohtsuka, G.; Glueck, J.A.; Mueller, J.; Nose, Y. Hemolysis test of a centrifugal pump in a pulsatile mode: The effect of pulse rate and RPM variance. Artif. Organs 1997, 21, 1284–1287.

46. Tayama, E.; Nakazawa, T.; Takami, Y.; Makinouchi, K.; Ohtsubo, S.; Ohashi, Y.; Andrade, A.J.; Glueck, J.; Mueller, J.; Nose, Y. The hemolysis test of the Gyro C1E3 pump in pulsatile mode. Artif. Organs 1997, 21, 675–679.

47. Qian, K.X. Pulsatile impeller heart: A viable alternative to a problematic diaphragm heart. Med. Eng. Phys. 1996, 18, 57–66.

48. Wang, S.S.; Chu, S.H.; Chou, N.K.; Qian, K.X. The pulsatile impeller pump for left ventricular assist. Artif. Organs 1996, 20, 1310–1313.

49. Farrar, D.J.; Bourque, K.; Dague, C.P.; Cotter, C.J.; Poirier, V.L. Design features, developmental status, and experimental results with the Heartmate III centrifugal left ventricular assist system with a magnetically levitated rotor. ASAIO J. 2007, 53, 310–315.