Nonlinear Closed-Loop Predictive Control of Heart Rate and Blood Pressure using Vagus Nerve Stimulation: An In Silico Study

A. Cardiovascular Model

The rat cardiovascular system (CVS) model, adapted from [18], includes five compartments representing the left heart, the arteries and the veins in the upper and lower body. The right heart and the pulmonary circulation is modeled by a capacitance, which is added to the venous capacitance in the upper body.

The model mimics an electrical RC-circuit (see Fig. 1). The blood pressure in the four arterial and venous compartments, and the total volume of the left heart follow Kirchhoff’s law:

where the abbreviations of the subscripts are listed in Table I. $R$ is resistance; $P$ is pressure; $V$ is volume; $C$ is constant compartment compliance. Pumping is achieved by a time-varying elastance of the left heart, which transitions between an exponential function representing diastole and a linear function representing the end of the systole [19]:

where $E_{max,lh}$ is the maximal contraction elastance of the left heart; $V_{u,lh}$ is the unstressed volume of the left heart; $P_{0,lh}$ and $k_{E,lh}$ are constant parameters for the exponential function. The term $\varphi (t)$ represents the periodic “activation function”. When $\varphi (t)=1$, the left heart is at maximum contraction, when $\varphi (t)=0$, the left heart is at complete relaxation. Finally, the heart valves are modeled as diodes, which results in four different time intervals in a normal cardiac cycle.

B. Baroreflex Model

The baroreflex model is adapted from [19]–[21]. It consists of the baroreceptor, afferent pathway, efferent pathway, and the effectors in CVS (see Fig.2).

The baroreceptors are stretch-sensitive fibers located in the carotid sinus which convert the arterial pressure to afferent firing rates. A model with two Voigt bodies and a spring in series is used to describe the stretching dynamics of the baroreceptive nerve endings:

where ${ϵ}_1$ and ${ϵ}_2$ are the relative displacements within each Voigt body; ${\alpha }_1$, ${\alpha }_2$, ${\beta }_1$, and ${\beta }_2$ are nerve ending constants; ${ϵ}_{ne}$ denotes nerve ending strain; ${ϵ}_w$ refers to the strain of the arterial wall which is described by a static sigmoid function of arterial pressure. The firing rate in the afferent pathway is calculated by a leaky integrate-and-fire model represented by:

where $I_{ne}$ is the current stimulus injected into afferent fibers, modeled as a linear function of ${ϵ}_{ne}$; $g_{leak}$ is a leakage conductance; $C_m$ denotes the membrane capacitance; $V_{th}$ is a given voltage threshold, and $t_{ref}$ is a refractory period.

The afferent firing rate is translated into (1) an efferent sympathetic firing rate using a monotonically decreasing function; and (2) an efferent vagal firing rate using a monotonically increasing function with an upper saturation, respectively. These efferent signals are delivered into the end organs in the CVS. The regulation effects include changes in peripheral resistance in upper and lower body, heart contractility, total stressed blood volume and heart period. The response of the heart period includes a balance between the sympathetic and vagal activities. The heart period changed by sympathetic stimulation ($T_{es}$) includes a time delay, a logarithmic static function and a low-pass first-order dynamic, while the response of heart period to vagal activities ($T_{ev}$) shows a positive linear relationship. Finally, the heart period is achieved by linear interaction between sympathetic and vagal responses. The dynamics for changes in heart period are given in (11)–(13). Dynamics of peripheral resistance, heart contractility, and total stressed blood volume are similar to those of $T_{es}$.

where ${\tau }_{T_{es}}$, $G_{T_{es}}$ and $D_{T_{es}}$ (${\tau }_{T_{ev}}$, $G_{T_{ev}}$ and $D_{T_{ev}}$) are the corresponding time constant, gain and delay in sympathetic (vagal) pathway. $T_0$ is intrinsic heart period.

C. Device Model

A device model is developed to predict the response of firing rate of different fibers to VNS pulse trains. In the rat, the barofibers are bundled into the aortic depressor nerve, which follows alongside the vagal nerve within the vagal nerve bundle [4]. The cervical sympathetic nerve runs separately from the vagus in rats. However, other mammals, including dogs, pigs and humans have a cervical vagosympathetic trunk [22], [23] and detailed microdissection to identify and isolate the sympathetic trunk, the aortic depressor nerve, and the vagal nerve may not be safe or practical during human VNS surgery [24]. There is also a big difference in cardiovascular control between the left and right VNS. In rats, left-sided cervical VNS caused stronger bradycardiac, hypotensive, and tachypneic effects than right-sided VNS [25]. However, this contrasts with findings in dogs [26] that demonstrated greater bradycardia with right-sided than left-sided VNS.

The device model we propose predicts rat cardiovascular response to left VNS. It engages three fiber types during VNS, representing the barofibers, the vagal fibers, and the sympathetic fibers. Each type of fiber distributes nonhomogenously in different stimulation locations (Fig.3, Left) and the total firing rate of each type of fiber is calculated as the weighted sum of its firing rate in each stimulation location:

where $f_i$ is total firing rate of fiber type $i$; $n$ is number of stimulation locations; $L_j$ represents the on-off condition of location $j$; $C_{ij}$ is concentration of fiber type $i$ in location $j$; $f_{ij}$ is firing rate of fiber type $i$ in location $j$; $f_{i,phy}$ is physiological firing rate of fiber type $i$ unaffected by external stimulation.

The proposed device model for activating each fiber type in each stimulation location (Fig.3, Right) consists of a nerve fiber activation curve and a conduction map. The model assumes that the width of the stimulation pulse train is constant, and that the pulse amplitude and frequency varies. The activation curve represents the proportion of nerve fibers recruited by stimulation of a particular amplitude using a static sigmoid function [27].

where $P_{ij}$ refers to activation probability of fiber $i$ in stimulation location $j$; $I_{j,stm}$ is the amplitude of VNS pulse train; $I_{ij,mid}$ is the value of amplitude at central point of the sigmoid function, and $k_{I_{ij}}$ is the slope at the central point. A conduction reliability (${\lambda }_{ij}$) represents the percentage of total action potentials conducted to the somatic end of the nerve fiber by considering the interactions between stimulation and physiologically induced action potentials. As the frequency of the stimulus and/or rate of the physiological input increases, conduction reliability decreases due to collisions between these two types of action potentials, as well as the inter- and intra-signal loss of excitability due to neural refractory period. We adopted the conduction map from [28] by fitting the data from the mechanistic model to a linear combination of physiological frequency and external stimulation frequency. The relationship between stimulation parameters and fiber firing rate at each stimulation location can be represented by the following equation:

where $f_{j,stm}$ refers to firing rate of the external stimulus at location $j$; ${\lambda }_{ij}^w=P_{ij}{\lambda }_{ij}$ refers to weighted conduction reliability.

D. Model Parameter Estimation

Since the model was based upon several compartment models previously developed in [18]–[20], we either used the same parameter values as those used in the original models or adapted some values to produce simulations that best matched experimental data associated with nominal, diseased, and exercise states of rats. We describe below the various steps used in this model validation. The results of the model validation and matching with experimental data are provided in section IV A.

To simulate the nominal states, the parameters in the CVS model were adapted from those in [18] to match the experimental data of young healthy rats in [29]. The parameters in the baroreceptor model are the same as those used in [20]. The parameters in the efferent pathway are the same as those used in [19], while the parameters in the transfer function related to effector dynamics are modified by multiplying the original value by a scaling constant which indicates the difference in heart period and blood volume between a rat and a human. Three stimulation locations are identified for the device model. The parameters representing the fiber concentration, the activation curve and the conduction map in baroreceptive and non-baroreceptive locations are obtained by matching the simulated change of HR and MAP by VNS to experimental data in [4], [30], respectively. Percent of vagal fibers which are not contained in these two locations are set to the third stimulation location.

Hypertensive heart disease is used as an example to show the feasibility of closed-loop control because VNS has an acute effect on reducing hypertension. Other cardiovascular diseases responsive to VNS require activation of long-term mechanisms which are not captured by the model. Drug-resistant hypertension is related to increased arterial stiffness, vascular remodeling, and increased sympathetic activity and decreased vagal activity. Most longstanding hypertension ultimately leads to heart failure and diastolic dysfunction is characterized as an early marker of heart damage in hypertension. To simulate the diseased state, we introduce an offset in sympathetic and vagal activity, increase the cross-section area of the artery in relaxation state and modify the gain of each effectors to match the data in [31].

Acute exercise triggers multiple physiological responses, including redistribution of blood flow and modification of sympathetic and vagal activities by central command. To address these issues, we separate the peripheral resistance in the lower body into two parallel resistors: $R_{sl1}$ representing the resting muscle resistance and $R_{sl2}$ representing the active muscle resistance. The exercise states are determined by a combination of an offset to sympathetic activity and forced change in active muscle resistance ($R_{sl2}$). The experimental hemodynamic data of Sprague-Dawley rat during last minute of treadmill exercise from [32] is used as reference for tuning the parameters for modeling the nominal rat. When estimating parameters for an exercised hypertensive rat, we also selected different values for gain parameters of the effectors to match the arterial pressure and HR experimental data for spontaneous hypertensive rat (SHR) during the last minute of the treadmill exercise in [33].

A. Reduced Model

The reduced model has the same components as the previously described full simulation model. The primary difference is that the reduced model only considers the inter-beat dynamics and ignores the intra-beat dynamics. Some simplifications of the cardiovascular and baroreflex model are made, while the device model is kept the same as the full simulation model. The CVS model is reduced to a three-compartment model, representing the artery, the vein, and the left heart. The time-varying elastance of the pumping heart is replaced by the cardiac output as a function of heart period, venous and arterial pressure using the integration method developed in [18]. The reduced cardiovascular model can be finally represented by the following differential equations:

where $P_{vc}$, $P_{ao}$ are the pressure of the veins and the arteries, respectively; $C_{vc}$, $C_{ao}$ are the compliance of the veins and the arteries, respectively; $R_{sys}$ is the systemic resistance in the peripheral circulation; $Q$ is the cardiac output; $E_f$ and $E_e$ are the averaged elastance of the left heart during filling and ejection, respectively.

The previously described baroreflex model is sensitive to the functional form of arterial pressure and includes a time delay. But the arterial pressure predicted by the reduced cardiovascular model is constant during each cardiac cycle. As a result, a reduced baroreflex model is adapted from [34] to regulate the averaged arterial pressure. The physiological firing rate of the barofibers ($f_{as,phy}$) is defined as the product of the ratio of the arterial pressure to some predefined setpoint of the arterial pressure ($P_{sp}$) and a predefined setpoint of the baroreceptive firing rate ($f_{as,sp}$).

The firing rates of sympathetic fibers ($f_{es,phy}$) and vagal fibers ($f_{ev,phy}$) are related to the baroreceptive firing rate by the following monotonically decreasing and increasing static curves, respectively.

where $f_{es,max}$, $f_{ev,max}$ are maximum firing rates of sympathetic and vagal fibers; $k_{es}$ and $k_{ev}$ are steepness parameters of sympathetic and vagal pathway; $\delta =F\left(f_{as,phy},u\right)/f_{as,sp}$ is the normalized baroreceptive activity modified by external stimulation with $F(\ast )$ representing the device model defined in (14) – (16) and $u={\left\{L_j,I_{j,stm},f_{j,stm}\right\}}^T$ is the vector of control variables representing the on-off conditions of three stimulation locations and the corresponding amplitude and frequency of the pulse train. Five effectors, representing the heart period, the systemic resistance, the elastance of left heart during filling and ejection, and the total unstressed blood volume are responsive to efferent drive. The change of each effector (${\theta }_i$) is described by the following [35]:

where $n_{es}=\frac{F\left(f_{es,phy},u\right)}{f_{es,max}}$ and $n_{ev}=\frac{F\left(f_{ev,phy},u\right)}{f_{ev,max}}$ are normalized firing activities of sympathetic and vagal fibers with external stimulus, respectively; ${\alpha }_i$, ${\beta }_i$, and ${\gamma }_i$ are the gain parameters; ${\tau }_{{\theta }_i}$ refers to the characteristic time constant. This uniform formulation removes the discontinuity in a logarithmic function and makes the nonlinear optimization tractable.

B. System Identification

Several of the parameters in the reduced model need to be recalculated to match the dynamics of the full in silico simulation model. The nonlinear grey-box technique was used to identify the reduced model with the MATLAB implementation of the algorithm ‘nlgreyest’. The selected parameters to be identified are listed in Table II; the other parameters are kept the same as in the full simulation model. The algorithm estimates the reduced model parameters by minimizing the error between the reduced model output and the output from the full model using the following objective function:

where $t$ is time, $N$ is the number of data samples, and $\epsilon (t,\eta )$ is the error between the reduced and full model output; $\lambda$ is a positive constant which trades-off variance versus bias error; $R$ is the weight matrix for variance error; ${\eta }_{min}$, ${\eta }_{max}$ are vectors of the lower and upper bound of each parameter which is selected to be ±5-fold of the corresponding parameter value in the full model. The input-output signals are generated from the full diseased model in both rest and exercise regimes. A uniform distribution is used to generate 50 random values of pulse amplitude and pulse frequency within their range for each of the 8 combinations of stimulation locations. Each perturbation of stimulation configurations lasts for 20 cardiac cycles to capture the dynamic response of the system. The output signals are the cycle-averaged mean arterial pressure (MAP) and heart rate (HR) and the sample time was set to 1 for integration.

The reduced model predictions of the diseased condition in rest and exercise regimes are shown in Fig. 5 as dashed lines, and the data generated by the full diseased model is shown as solid lines. The reduced model has a 85.97% match for the MAP and 92.54% match for the HR for rest regime (and 75.27% and 84.5% for exercise regime). The mismatch between the identified reduced model (“predictor”) and the full in silico model can be expected to be handled by the feedback action of the NMPC with the moving horizon estimator.

C. Nonlinear Model Predictive Controller

NMPC is based on an optimal control algorithm. At each sample time, it calculates several control actions over a future time horizon by minimizing an appropriate cost function over a future prediction horizon using the reduced model to predict the system response to these control moves. Only the first control action is applied to the system, new measurements of HR and MAP are obtained and the prediction and control horizons are shifted one step forward to compute the next set of optimal control moves.

The objective of NMPC in our VNS context is to bring HR and MAP to its nominal value in both rest and exercise scenarios. Previous studies on closed-loop VNS systems used either the stimulation frequency or amplitude as the only manipulated input, but all stimulation parameters, as well as stimulation locations, affect therapeutic efficacy, and should be considered simultaneously when making control decisions. The HR and MAP signals are sampled by the controller at the end of each cardiac cycle. The NMPC algorithm is designed to calculate control actions during the cardiac period and deliver the control decision to the CVS synchronized with each cardiac cycle.

The NMPC consists of three elements: the estimator, the regulator, and the target calculator. Each element solves an optimization problem by incorporating the reduced model with given constraints on input or state variables. The regulator formulation is as follows:

where $k$ is the current cycle index; $N$ is the prediction horizon; $x_s$ and $u_s$ are steady state value of the states and inputs, which are calculated by a standard target problem for each setpoint of the output; $\text{Δ}{u}_k$ is the change in input variables $u_k$; $Q_r$, $R_r$, $R_u$, and $P_r$ are weighting matrices. The first and second terms in the objective function penalize the deviations of the state predictions from their reference. The third term penalizes large input changes, and the last term provides a terminal constraint to achieve closed-loop stability, where $P_r$ is calculated using the Riccati equation by linearizing the system around the steady-state point. The first two equality constraints represent the reduced model as discussed above, in discrete form. $\widehat{x}$ is the future state predicted using the reduced model; $\widehat{d}$ is the estimated value of the unmeasured disturbance; $B_d\in {ℝ}^{N_u\times N_d}$, $B_{pd}\in {ℝ}^{N_p\times N_d}$, and $C_d\in {ℝ}^{N_y\times N_d}$ are input matrices for disturbance. The third constraint imposes a lower bound and upper bound on the inputs to ensure an appropriate stimulation intensity. The last constraint is a binary equality constraint for stimulation locations – each stimulation location is in closed state (value 0) or open state (value 1).

The resulting regulator problem can be treated as a mixed-integer nonlinear program (MINLP). Solving a MINLP is a computationally difficult task since the optimality and speed of convergence cannot be guaranteed. Instead of applying an MINLP solver, we explicitly list the 8 possible combinations of the 3 stimulation locations. For each feasible combination, the explicit values of the 3 binary input variables are inserted into the regulator formulation to give a traditional nonlinear program without binary variables. By repeating the calculation for each feasible combination, the one with the minimal cost is chosen as optimal.

In practical situations, physiological parameters may change with time, which causes unanticipated noise and disturbance on output variables. A moving horizon estimator (MHE) is developed to address these issues using a series of measurements observed over time to estimate the current state and disturbance. The estimation problem is stated as:

where $N_t$ is the length of the moving horizon window; $X={\left[x_0,x_1,\dots ,x_{Nt}\right]}^T$ represents the sequence of estimated state variables; $W={\left[w_0,w_1,\dots ,w_{Nt}\right]}^T$, $V={\left[v_0,v_1,\dots ,v_{Nt}\right]}^T$, and $D={\left[d_0,d_1,\dots ,d_{Nt}\right]}^T$ represent the sequence of estimated state noise, measurement noise, and disturbance, respectively; $P_e$, $P_{d_e}$, $Q_e$, $Q_{d_e}$, and $R_e$ are weighting matrices. The first and second terms in the objective function penalize the arriving cost of state and disturbance variables, respectively. The third and fourth terms minimize the process and measurement noise, respectively. The last term prevents large change in subsequent disturbance estimation. Similar to the regulator problem, the MHE is solved repeatedly at each sampling instance. At each estimation step, the first element from the previous estimation is eliminated and the newest measurement is added into the estimation window. The full simulation model does not introduce any process and measurement noise so that the $Q_e$ and $R_e$ are set to a very small value in the optimization problem. By estimating the unmeasured disturbance of model parameters, the steady-state calculation may detect an infeasible problem. In this case, the objective function of the regulator is modified to the following standard MPC formulation:

The computational cost for solving an optimization problem in NMPC depends on its size and formulation. While large horizons are often necessary to capture slow dynamics, the large horizon also significantly increases the computational cost. Thus, a trade-off between optimal solution quality and computational cost is considered in the NMPC design. Here, the control, prediction and the estimation horizons are all set to 10 cycles, which are large enough to allow the controller to provide satisfactory performance.

A. Open-Loop Model Validation with Data

The open-loop behavior of the full in silico simulation model is compared to experimental data from the literature.

B. Evaluation of NMPC

The NMPC algorithm is evaluated using three simulation scenarios: 1) tracking of nominal state, 2) disturbance rejection, 3) adaptation to exercise.