## Abstract

Inverse problems in cardiovascular modelling have become increasingly important to assess each patient individually. These problems entail estimation of patient-specific model parameters from uncertain measurements acquired in the clinic. In recent years, the method of data assimilation, especially the unscented Kalman filter, has gained popularity to address computational efficiency and uncertainty consideration in such problems. This work highlights and presents solutions to several challenges of this method pertinent to models of cardiovascular haemodynamics. These include methods to (i) avoid ill-conditioning of the covariance matrix, (ii) handle a variety of measurement types, (iii) include a variety of *prior* knowledge in the method, and (iv) incorporate measurements acquired at different heart rates, a common situation in the clinic where the patient state differs according to the clinical situation. Results are presented for two patient-specific cases of congenital heart disease. To illustrate and validate data assimilation with measurements at different heart rates, the results are presented on a synthetic dataset and on a patient-specific case with heart valve regurgitation. It is shown that the new method significantly improves the agreement between model predictions and measurements. The developed methods can be readily applied to other pathophysiologies and extended to dynamical systems which exhibit different responses under different sets of known parameters or different sets of inputs (such as forcing/excitation frequencies).

## 1. Introduction

Numerical models of blood flow of varying complexity—ranging from *lumped* models that are electric analogues of the circulation, through *one-dimensional* (1D) Euler equations and *three-dimensional* (3D) Navier–Stokes models, to *geometrical multi-scale* models—have been developed to study various healthy and diseased physiologies. The lumped parameter models, in particular, find application either individually or as boundary conditions (open or closed loop) to higher order models (1D or 3D) in the geometrical multi-scale method [1,2]. For a patient-specific analysis, however, these models must be adapted to each patient individually. This means that the model parameters must be estimated from some set of clinical measurements (typically noisy) in the patient. This article deals with important perspectives related to such an estimation when data are assimilated gradually during the simulation, i.e. in a sequential data-assimilation setting.

Estimating model parameters in haemodynamics is challenging due to three factors: the presence of uncertainty (noise or biases) in clinical measurements, the imperfect nature of various assumptions employed in the models and the fact that the number of parameters can be quite large even in relatively simplistic models representing the whole circulation. The inverse problem has been addressed in many previous studies, see for example [3–16] and [17,18] for their advantages and disadvantages. A recent Bayesian approach to parameter estimation in cardiovascular models is based on Markov chain Monte Carlo (MCMC) methods with mean or maximum target quantities [19]. In recent years, the data-assimilation approach has gained popularity as it can take advantage of full-time varying measurements as opposed to just scalar mean/minimum/maximum indices. It has been employed in a range of cardiovascular problems such as the estimation of Windkessel parameters [18,20–22], closed-loop lumped model parameters for single-ventricle physiology [17], artery wall stiffness parameters [22,23] and tissue support parameters for fluid–structure interaction models [24]. Recently, it has also been used for parameter estimation in electrocardiography/electromechanical models [25–28], cardiac motion recovery [29,30], inverse electrophysiology applications [31,32], and for baroreflex regulation modelling [33]. This study focuses on some aspects of such a data-assimilation method, namely the unscented Kalman filter (UKF), towards successful estimation and interpretation in closed-loop lumped parameter models of cardiovascular haemodynamics.

The first contribution of this study is the management of the Sigma point collapse in the sequential UKF framework which can lead to ill-conditioning of the empirical covariance matrix, and consequently numerical problems. The second contribution of this study is to present a method to include the measurements of minimum (min) and maximum (max) values of model outputs in the method in an efficient manner. This is relevant as sometimes only min/max values, as opposed to full time-varying measurements, are available/measurable for certain quantities in the clinic. The third contribution is the presentation of different forms of prior knowledge that can be incorporated in the UKF setting to ensure that a reasonable solution of the inverse problem is obtained. The final contribution is a method to include measurements at different heart rates, which is a challenge in the basic UKF method.

The first two contributions are presented in two patient-specific cases of congenital heart disease, namely stage I single-ventricle physiology (one with normal valves and one with a regurgitant atrioventricular valve). The last contribution of assimilating measurements at variable heart rates is demonstrated first for a synthetic dataset, to verify the approach, and then on a patient-specific case. This method can be readily extended to other dynamical systems that exhibit different responses under: (i) differing values of a single (or a set of) known parameters (i.e. parameters that are not to be estimated), for example estimation of medium diffusivity and conductivity in an advection–diffusion problem where species concentrations are measured at two different known advection velocities; and (ii) different inputs, for example different input signals or different forcing/exciting frequencies. The latter case is especially relevant as, in practice, certain models can be particularly sensitive to some parameters in certain time zones (or space–time zones in the case of partial differential equations) [34]. It is necessary to gain information either by having very precise information in these zones or by visiting them several times sequentially, so that they are revisited when knowledge on the other parameters has been gained. Assimilating data over several quasi-periodic signals is such an example that is naturally encountered in physiology.

It should be noted that single-ventricle physiology presents a challenging treatment as only one heart pump functions well at birth; treatment is usually performed in three stages of open-heart surgery where important decisions must be made at each stage [35,36]. In recent years, computational fluid dynamics (CFD) models, large parts of which comprise lumped parameter models, has been employed to aid clinical decision-making and development of novel surgical strategies [36,37]. For individual case-by-case clinical decision-making in single-ventricle patients, however, the estimation of lumped model parameters while taking into account the uncertainty in clinical measurements is necessary. It is this problem that the methods developed in this article relate to.

The model of single-ventricle physiology adopted in this study is presented in [17] and employs the one-fibre model [38,39] to describe the heart chambers and an extension to the valve model of [40] presented in [17]. The main motivation for such model choices is that they are based on physically meaningful parameters such as wall volumes of the heart chambers, effective areas of the valves, sarcomere properties, etc., for which clinical/experimental estimates are available. This reduces the black-box nature of the lumped models, and routine measurements such as four-chamber electrocardiography, Doppler velocimetry, etc., can be employed, at least qualitatively, for validation of some of the parameter estimates. This moreover allows the estimated lumped model parameters of the heart to serve as a first good guess for myocardium material properties in larger 3D models.

## 2. Material and methods

### 2.1. Single-ventricle circulation model

The model for the entire circulation (so-called closed-loop model) for the single-ventricle case after stage I surgery considered in this study is described in [17]. Here it is briefly presented. An electric network schematic of this physiology along with the pressure–volume relationships of various model components is shown in figure 1. Here, the resistances ( for linear dependence of pressure on flow rate and for quadratic dependence) represent the viscous losses in the blood; the capacitances ( ) represent the compliance of arteries or veins due to the elastic nature of their walls; and the inductances ( ) represent the fluid inertial effects. Note that an artificial shunt connects the systemic and pulmonary circulations at this stage of the heart palliative procedure. The models for the heart chambers ( ) and the valves, which function to separate the heart chambers from the rest of the circulation, ( ), are described next; these, combined with mass conservation at each network junction of figure 1, result in a system of nonlinear ordinary differential equations (ODEs).

#### 2.1.1. Heart chambers

The heart chambers are described by a one-fibre model [38,39], where the relationship between wall fibre stress *σ*_{f}, pressure inside the cavity *p* and volume of the chamber *V* is described by
2.1where *V*_{w} is the wall volume of the chamber. The fibre stress is the sum of an active component *σ*_{a} responsible for the heart contraction during each heart beat and a passive component *σ*_{p}. Let *V*_{0} and *l*_{0} represent the volume of the chamber and the length of the sarcomere, the functional unit of the fibre, respectively, at zero transmural pressure when there is no contraction. The fibre stretch *λ* relating the sarcomere length *l* to the volume *V* of the chamber is given by
2.2 Consequently, the sarcomere shortening velocity, , can be computed as . The active component of the stress is described as
2.3
2.4
2.5 and
2.6where , , and describe the dependence of sarcomere length *l* on active stress [17], is the time elapsed since contraction activation, is the activation duration of the chamber in one cardiac cycle, is the initial sarcomere shortening velocity, is the maximum active sarcomere stress and is a shape parameter. The passive stress is given by
2.7 where and are sarcomere material properties.

Similar to [17] is assumed and is parametrized by three parameters: , the activation duration of the heart chambers, namely the single atrium; , the activation duration of the single ventricle; and , the overlap between the end of atrium activation and the beginning of ventricle activation (figure 2).

#### 2.1.2. Valve description

The pressure drop across a valve is given by the following Bernoulli relation:
2.8where *ρ* is the density of blood, and and denote the effective opening area and effective length of the valve, respectively. The valve state is described by a single variable, , that relates to the effective area at time *t* as follows:
2.9where and is negative only if the valve is regurgitant due to prolapse, and the parameters , and denote the maximum, minimum and maximum regurgitant area under prolapse, respectively. The valve dynamics is described by ODEs such that the rate of valve opening/closing is proportional to the favourable pressure difference (see [17] for details)
2.10where is the pressure gradient beyond which the valve prolapses, and are proportionality parameters.

### 2.2. Overview of the unscented Kalman filter

After having described the hemodynamics model above, let us explain how its parameters are estimated based on dynamics measurements of blood flow and pressure. The UKF is an extension to the Kalman filter for nonlinear problems [41,42]. Let be the state vector. The so-called *forward* model that relates the state at time to that at can be written as
2.11where represents the parameter vector, represents the time step, and *F* represents the forward operator. For example, in the model represented in figure 1, the state vector is represented as . The operator *F* is then defined by the discretization of the ODE system resulting from the heart and valve models described in §§2.1.1 and 2.1.2 and the individual component equations shown in figure 1 (right). The measurements are described by the observation operator as
2.12where is the observation vector, i.e. the subset of measurements chosen to estimate parameters, at time , *H* represents the observation operator that relates the model state to the measurements, and represents the measurement noise. For parameter estimation through sequential filtering methods, the state vector is combined with the parameter vector to yield an augmented state vector , and the forward, , and observation, , operators are written as
2.13and
2.14 The goal of sequential methods is to provide an estimate for the state by taking all the measurements to into account. This is achieved recursively by providing an initial estimate with a covariance matrix . The recursion is as follows: first, the state and covariance are propagated to an intermediate state with covariance through the forward model of equation (2.13), and then the estimate with covariance is obtained by correcting the intermediate state, by accounting for (assimilating or probabilistically conditioning on) the observation through equation (2.14). If () observations are available, then an estimate of the parameters is provided by the corresponding components of , and associated variances by the corresponding diagonal elements of .

If the forward and observation operators for the augmented state and are linear, and is white noise, then the Kalman filter [43] provides an optimal state estimate and both the propagation and assimilation steps can be analytically computed. For non-linear problems, various extensions such as the extended Kalman filter (EKF), ensemble Kalman filter (EnKF), UKF, etc. have been proposed, which differ primarily in the forward propagation step. For a discussion on such choices, the reader is referred to [18]. In what follows the forward propagation and assimilation steps of the UKF are briefly described.

#### 2.2.1. Forward propagation of uncertainty

In UKF, forward propagation is performed by propagating a set of deterministically chosen *particles*, called *Sigma points*, through the forward model of equation (2.13). The Sigma points, denoted by , (), are generated as follows:
2.15where and represents the *i*th row of the matrix square root of . Each is then propagated as
2.16where denotes the propagated particle. and are calculated from the propagated Sigma points as
2.17and
2.18where and represent the weights associated with Sigma points for empirically calculating the mean and covariance, respectively [18]. In a similar manner, the intermediate (prior) mean and covariance of the observation vector are calculated by propagating through equation (2.14) and calculating the empirical statistics
2.19
2.20
2.21where represents the covariance matrix of the measurement noise . Finally, the cross-covariance matrix between and is
2.22

#### 2.2.2. Assimmilating the observations

Through the measurement , the above-calculated intermediate state is corrected to yield an estimate of the state with covariance at by minimizing the following cost function : 2.23

The solution to the above in the UKF setting is given by 2.24and 2.25where the gain matrix is given by 2.26

### Remark 2.1.

It should be noted that if the observation operator is linear then the propagation of the Sigma points through the observation operator, equation (2.19), is not necessary. Instead, and in equations (2.21) and (2.22), respectively, can be analytically computed. For ODE models, in certain cases, this may induce a choice of formulation. Consider the heart model as an example. This model (see figure 1 and §2.1.1) is highly nonlinear where, for each heart chamber, the ODE is formed naturally for volume *V* as follows:
2.27and the pressure, *p*, is calculated by equations (2.1)–(2.7). However, it is the pressure, or , that is typically observed. If one were to reformulate the ODE in terms of pressure (by differentiating equations (2.1) and (2.27)) rather than the volume, then pressure would be a part of the state and hence a linear observation manager would suffice. However, given the nature of equations (2.1)–(2.7), it is much easier to keep the ODE state variable as *V* and employ a nonlinear observation operator.

### 2.3. Open loop versus closed loop

In the context of estimating parameters for haemodynamic systems, such as the one presented in figure 1, a natural question that arises is whether one should perform parameter estimation in the entire closed loop or perform simpler estimation in sub-segments of the closed-loop circulation, i.e. in an open-loop setting. In terms of computational cost, solving multiple smaller open-loop segments is cheaper than one larger closed-loop model, and this may be a driving motive for open-loop estimation. However, in ODE systems such as that of figure 1, the cost of solving the full closed-loop model is not exorbitantly high, and hence the computational gain would be of little practical significance. Another consequence of estimation in open-loop segments is related to the manner in which uncertainty in the measurements is accounted for. As an example, consider the lower body segment in figure 1 (the region between nodes N-1 and N-2). As at both ends, the flow-rates, i.e. and , and the pressures, i.e. and , are available as measurements, one may consider an open-loop model just for this segment. However, in order to solve the forward model one must impose two forcing conditions, in this case say at the ends N-1 and N-2. For this, one may choose either *p* or *q* at each end and employ the remaining two measurements as observations to the UKF. While this is a reasonable approach, in doing so two of the measurements at the boundaries have been deterministically imposed and hence are considered devoid of any uncertainty. This might lead to erroneous results if such an assumption, depending on the quality of the measurements, cannot be justified. Furthermore, there is little assurance that, when the estimated parameters for each segment are plugged into the closed-loop model with no such forcing conditions, the model output will agree with the measurements. On the other hand, there is a distinct advantage in an open-loop setting. In certain cases, the measurements are not synchronized in time with respect to each other. This is especially true if simultaneous electrocardiogram (ECG) measurements are not taken or are unreliable, or due to different techniques for the pressure (typically catheterization) and flow-rate (typically magnetic resonance imaging (MRI)) measurements, implying that their simultaneous measurement is rarely possible. In such cases, some of the time shifts can be posed as parameters in the open-loop models, and can be estimated by the UKF. As an example, consider the aforementioned lower body open-loop segment. Assume that the pressure measurement is synchronized with and the flow-rate measurement is synchronized with , but the pressure and flow rate are not synchronized with respect to each other. In the UKF, one may impose the following as forcing conditions (FCs) and measurements while accounting for the phase/time shift between and by a parameter
2.28and
2.29Using the above formulation, and by including as one of the parameters to be estimated, the discrepancy due to asynchronous measurements can be taken into account. Note that the parameters must affect the forward model in order to be estimated; hence only in an open-loop setting where appears in an FC and this affects the forward model can it be estimated. Such an application of estimating the asynchronization is presented in [18]. By contrast, in a closed-loop setting the asynchronization between the measurements cannot be automatically accounted for.

Preliminary tests suggested that parameter estimation on individual open-loop segments of figure 1 may not result in as good a reproduction of all the measurements when substituted in the closed-loop model (also see [44] where in an open-loop configuration the heart parameters estimated for patient A are considerably different from those estimated in a closed-loop setting in [17]). Thus, in what follows, estimation is performed directly on the closed-loop model. Such an observation suggests that uncertainty in the measurements is an important consideration.

### 2.4. Sigma point collapse

The valve model presented in §2.1.2 is essentially described by the valve dynamics parameter *ξ*. In haemodynamics problems, the valve state remains closed or open for long periods of time in a cardiac cycle. This implies that during these periods the valve state *ξ* remains 0, 1 or −1, and all perturbations (i.e. Sigma points) to this variable at a time may result in the same final state (0, 1 or ) at time . In two dimensions, this phenomenon is depicted in figure 3, where at all Sigma points have the same value for the *ξ* variable. This phenomenon makes the empirical covariance matrix not positive definite, which may result in a non-positive definite after the data-assimilation step is performed. Note that the matrix square root in equation (2.15) is typically calculated via the Cholesky factorization, and hence the non-positive definiteness of presents problems. This issue can be remedied by adding a small diagonal term to in equation (2.18) as follows:
2.30where is an identity matrix and is a scaling factor. The added matrix can be interpreted as a regularization term, or as the covariance matrix of error in the forward propagation of the model. In the latter view, the discretized model can be seen, in lieu of equation (2.13), as
2.31where is the discretized model propagation error with zero mean and covariance . The end result is that, even though the Sigma points collapse in one or more dimensions, the propagated variance in these directions is maintained at a low value of . In the Bayesian perspective of UKF, the prior on can be seen as a normal distribution with mean and covariance . Then, the data assimilation step consists in computing the posterior distribution of . When collapse of Sigma points occurs, the prior on the collapsed states becomes a dirac delta function, i.e. the prior has no uncertainty associated with these states, and hence cannot be corrected in the posterior. The term avoids such a case by ensuring that no state variables become entirely certain at any step of the UKF.

### 2.5. Minimum and maximum of a time-varying output

While the UKF can easily incorporate time-varying measurements by time discretization, i.e. a measurement is provided at discrete time instants as observations, it is not straightforward to include measurements that correspond to the minimum or maximum of a time-varying output. This is a common measurement in haemodynamics problems, when only min/max volumes, pressures or flow rates are known. If the time instants when the minima or maxima occur are known *a priori*, then inclusion of such measurements is straightforward by only providing these observations at the known time instants. However, the locations of the maxima or minima depend on the parameters and the state variables, which are continuously being modified by the UKF. Therefore, methods to include such measurements are needed. One way, in a nonlinear observation operator setting, is to consider the observation operator that returns the max or min of an output, say the *i*th component of the state, as follows:
2.32where represents the *i*th component of , is the time step, *T* is the time period (one cardiac cycle), and
2.33

In the above method, a forward model needs to be solved for all the Sigma points for at least one cardiac cycle to obtain the min/max values of the observable output. In standard UKF, there are Sigma points, and consequently forward simulations are required. If the above method is utilized to include the min/max observations then, at each time when observations are assimilated, additional forward models need to be evaluated for at least one cardiac cycle. This significantly increases the computational cost associated with the UKF. As an example, if the UKF is run for two cardiac cycles and observations are assimilated at a total of *N* discrete time instants, then forward simulations (for two cardiac cycles) need to be run for the standard UKF, and additional forward simulations (for one cardiac cycle) to include the min/max observations. Alternatively, in certain cases, the time instants where the maxima/minima () for a particular output occur can be predetermined and are usually a function of the parameters. In such cases, the observation operator can simply be appended whenever to include the min/max measurement. As an example, consider that, for the ventricle, only the end-systolic volume (ESV) and end-diastolic volume (EDV) measurements are available. Typically, these will correspond to the maximum and minimum chamber volumes during a cardiac cycle. In the model, the contraction of the ventricle is dictated by the activation function parametrized by the parameters . Hence, given the parameter , the time instants where the ventricle will be at min/max volumes can be determined. This is shown in figure 4: the EDV (max) occurs just before the ventricle begins to contract, i.e. at , where is the time elapsed since the beginning of the current cardiac cycle, and the ESV occurs when the contraction phase of the ventricle ends, i.e. at . Hence, if the EDV and ESV measurements are available then one may consider a small interval of width around and , and provide the EDV and ESV measurements as the observations of only when the UKF time *t* is in these intervals, the intervals in turn being determined by the current state (which determines .

### 2.6. Prior knowledge: dynamical system and parameters

As the number of parameters to be estimated increases, the inverse problem of estimating parameters becomes harder. This is particularly true if the associated cost function between the model and the measurements, i.e. , is multi-modal or non-convex (i.e. it has either multiple minima or no minimum at all). In a classical least-squares minimization framework, identifiability is partly improved by adding a regularization term as follows: .2.34 As is evident from the above notation, the regularization term can be seen as the prior, i.e. the initial estimate for the state with mean and covariance , employed in the UKF method. As demonstrated in [18], a bad choice of and can lead to failure of the UKF method even in relatively simple dynamical systems. This aspect becomes more important when the size of the dynamical system and the number of parameters increase. Consequently, for such systems, prior knowledge on the parameters is important for the success of the method. The alternative is to run many instances of the UKF method with different priors, but this is a brute force approach that, although realizable, requires significantly more computational effort. Prior information may be classified into the following categories.

#### 2.6.1. Parameter values

Certain parameter values are known to be close to values determined clinically, through previous studies, or through experiments. For example, estimates of the wall volume of the ventricles, maximum valvular areas and activation durations of the heart chambers are available through clinical imaging, echocardiography or past clinical studies. In essence, this knowledge is a measure of what is expected in a population of patients of similar age and physiology, and is especially helpful when model parameters have a physiological meaning, such as effective valve areas, wall volumes of the heart chambers, etc. Such knowledge can be easily integrated in the UKF initial state . *A priori* there is generally no reason to assume correlation between the elements of and hence can be assumed to be diagonal with respective entries corresponding to the respective variances (inverse of confidences arising from the population-based knowledge) in the chosen prior estimates.

#### 2.6.2. Parameter ranges

A common form of prior knowledge is about parameter values in terms of known ranges of variation. Even more common is the knowledge that a certain parameter cannot be negative. For example all the resistances, capacitances and inductances in figure 1 are known *a priori* to be positive. Similarly, the volumes cannot be negative and activation durations cannot be larger than the cardiac cycle. While such information cannot be included in the form of a prior, one can include appropriate parametrizations to enforce such constraints. For positive parameters,
2.35 and for a parameter constrained between and
2.36may be chosen, where *Ψ* is the real parameter and the *ψ* is the form employed in the state . The parametrization ensures that the variable *ψ* can be manipulated by the UKF to take values from to in the Sigma points while yielding reasonable values of the real parameter to be employed in the forward model. The parametrization of equation (2.35) can also be employed for state variables such as pressure if they are required to be positive. This also makes sense from a conceptual viewpoint in the sense that, in the UKF, the measurement errors (variances diagonal elements of are seen as Gaussian, and such an assumption leaves a finite non-zero probability that the pressure might be negative. Instead, the parametrization of equation (2.35) transfers the normal assumption to the exponent *ψ*, thus making the real pressure *Ψ* log-normally distributed, which is a much more reasonable assumption due to only positive support of the log-normal distribution.

In the models used in this study, all pressures and volumes are parametrized by equation (2.35) to enforce positivity constraints. The state variables and are parametrized by equation (2.35) to be constrained between 0 and 1 or between −1 and 1 depending on whether regurgitation is included. Similarly, the parameters and are constrained between 0 and *T* (the time period of the cardiac cycle) and the overlap parameter is constrained between 0 and 0.05 s through equation (2.36).

#### 2.6.3. Prior knowledge on the range of variation of a state variable

The case when min/max values of certain state variables are known *a priori* with a reasonable degree of confidence is presented in §2.5. Here, the case when such ranges are not precisely known in a particular patient and yet there is some belief about what this range should not be is presented. For example, in stage I single-ventricle physiology, the ventricle volumes may vary from 5–25 ml during end-systole (end of the heart contraction period) to 15–50 ml during end-diastole (end of the cardiac chamber filling), depending on the age and physiology of the patient. A stronger belief is that these should not vary between 40 ml at end-diastole to 70 ml at end-systole, as these are too high for a three to six month old patient. This information may be seen as an imprecise knowledge about the mean value of ventricle volume. From equations (2.1) and (2.2), it is evident that, in the heart model, the ratios and determine the pressure. If only a noisy measurement of pressure in the aorta (reflective of pressure in the ventricle) is available, the parameters and are weakly identifiable: i.e. for a higher and a correspondingly higher , the aforementioned mean value of *V* might increase while yielding a similar value of pressure. The UKF, in such a case, may show a constant drift in both the parameters and and in the volume *V*. A natural question then is how to constrain the mean value from drifting too far from a known value, say . One way to include such a constraint is to provide the following measurement for the volume :
2.37where is a relatively high value in comparison with the real measurement variances. Essentially, for a time-varying quantity *V* this provides a constant pseudo-measurement of at each time of data assimilation. It results in a pulling force towards the value proportional to the distance of *V* from , and if appropriately chosen affects the solution only when drifts too far from . For example, assume that the true varies between 8 ml and 32 ml. If one provides a constant pseudo-measurement of 20 ml with a variance 400 ml (i.e. 0–40 ml with one standard deviation), the effect on the solution is negligible when the volume is in the true range, and will pull the solution towards 20 ml significantly when drifts far from 20 ml.

#### 2.6.4. Initial value of the state variables

A high degree of confidence can be associated with the initial values of some of the state variables based on the model physics. For example, in the heart model note that the driving functions are the activation functions shown in figure 2, and hence corresponds to the time just before the activation of the ventricle begins. Furthermore, *a priori* one may assume that the parameter is zero, implying no overlap between the atrium and ventricle activations, and hence at the atrium activation has just finished. Under these assumptions, confidence is high that at both the valves are closed, consequently both the valvular flows are zero, the aortic pressure is at end-diastole, the atrium is at end-systole and the ventricle is at end-diastole.

In a sense, all the above points are about choosing appropriate regularization for the inverse problem and how to include this regularization in the UKF method. These priors/regularizations are essential when estimating a high number of parameters and make the inverse problem better posed. Without appropriate prior information, one may see and the diagonal elements of as parameters of the UKF method. These are parameters, and searching for the right parameters to utilize in a method to estimate *p* parameters of the inverse problem defeats the purpose.

### 2.7. Heart rate variability in measurements

A common problem associated with haemodynamics measurements is heart rate variability. In particular, due to different techniques of measurements—and consequently different patient states during measurements (sometimes a difference of days)—the pressure (obtained by catheterization) and the flow-rate (obtained by MRI) measurements are typically at different heart rates. This poses problems in the basic UKF setting as the forward model is run only at one heart rate. A method to cope with such heart rate variability in the UKF is presented here. Consider the following two models at different heart rates: 2.38where and represent the state variables at time for the models and with heart rates and , respectively. Assume that the pressures are observed at and the flow rates at as follows: 2.39where are the observation operators, are the observations and are the corresponding measurement errors. For parameter estimation, the following augmented state models for forward and observation are considered: 2.40and 2.41Denoting the RHS of equation (2.41) by , the UKF is performed on the state with forward operator of equation (2.40) and observation operator of equation (2.41).

In the above formulation, the parameters are corrected by both the observations at HR_{1} and , but each state and is only individually corrected by the observations and , respectively. As most ODE systems for haemodynamics models are initial state forgetting, in the sense that irrespective of the initial state all runs eventually converge to the same state depending only on the parameters, it should be possible to estimate complete states and by correcting the parameters through observations at both heart rates. It should be noted that the above formulation requires forward and observation model simulations as opposed to in the basic UKF. Finally, the proposed method can be easily extended to cases where measurements at more than two different heart rates are available in a straightforward manner, albeit at a linearly increasing computational cost.

All the computations in this study are implemented in-house in the Python/Cython programming language. All simulations are run on an HPC node (RAM 48 GB) with two Intel Xeon X5650 processors, each with six cores and 12 threads.

## 3. Results and discussion

The results are divided into two categories: (i) when only a single heart rate is considered and (ii) when measurements are available at differing heart rates. For quick reference, a nomenclature of all the estimated parameters is presented in table 1.

### 3.1. Inverse problems considering a single heart rate

The model and methods for parameter estimation are applied to two patients: patient A with a normal atrioventricular valve (AVV) and patient B with a regurgitant AVV valve due to prolapse. Details of clinical data acquisition can be found in appendix A. Both patients correspond to stage I surgery for the hypoplastic left/right heart syndrome, a form of single-ventricle physiology. The parameter estimation results for these two patients are presented in [17] and the reader is referred to this study for the choice of model, time evolution of parameter estimates with the UKF method, validation of the parameter estimates and the clinical significance of the results. Here, the method to set up the UKF and the results of parameter variances, neither of which is included in [17], are presented. Parameter variances in particular can provide an indication of the identifiability of the system.

#### 3.1.1. Typical patient-specific case

For patient A, a nonlinear observation manager (see remark 2.1) is employed, Sigma point collapse is prevented for the valve states (§2.4) and the prior knowledge is considered as described in §2.6. In particular, in equation (2.31) is assumed to be with equal to . Additionally, to constrain the atrium and ventricle volumes, pseudo-measurements of and are respectively provided with ; see equation (2.37). Figure 5 shows the evolution of in the UKF for the 33 parameters. In this figure, which complements fig. 3 of [17], is the initial prior variance (diagonal elements of ) and is the UKF estimated variance (diagonal elements of ). Note that a negative value of in this plot corresponds to an uncertainty reduction by a factor of relative to the prior variance. The rate of decrease corresponds to the amount of information obtained from the measurements at each time: in general, a faster decreasing variance corresponds to a relatively easily identifiable parameter. Figure 5 shows that within four cardiac cycles, , the variance of all the parameter estimates, except and , is decreased by at least a factor of four relative to the prior variance. This trend of consistently decreasing variance continues as more measurements are added. By contrast, the parameters and appear to be weakly identifiable in the current setting as evidenced by the fact that in 20 cardiac cycles their variance has been less than halved. This may be related to small sensitivities of the measured model outputs to the parameters around the estimated values (see table 3 of [17] where the estimated values for these parameters are close to negligible). Small sensitivities and close to zero estimates indicate that the parameters are not critical to model output. By contrast, if they are critical to model output and appear to be weakly identifiable during estimation, it indicates a need for more measurements to improve identifiability. For example, the identifiability in the pulmonary block in figure 1 may be improved by providing an added measurement of the pulmonary artery pressure, . It should be noted that the idea of first identifying most influential parameters to the model output (or the measurements) and then proceeding with the estimation procedure has also been advocated in the literature; see for example [45–47].

Lastly, the parameter estimates obtained by the proposed method are clinically reasonable: for example, the ventricular wall volume measured by MRI for this patient is 34.2 ml, whereas that estimated by the method is 33.9 ml [17]. Similarly, the single-ventricle volume variation was measured via MRI between 10 ml and 29 ml, whereas that reproduced by the model is between 10.7 ml and 24.5 ml [17]. Measurements such as venrticular pressure which are not employed for parameter estimation also show a good agreement with the model output [17]. The parameter estimates and model output are further validated through other qualitative measures, such as echocardiography, Doppler velocimetry, *E*/*A* ratio for ventricular filling, and pulmonary venous wedge pressures; for details see [17].

#### 3.1.2. Enforcing ventricular volume extrema to model a case of valve regurgitation

Patient B includes atrioventricular valve regurgitation due to prolapse. In order to correctly estimate heart model parameters with a regurgitant atrioventricular valve, measurements pertaining to heart chamber volumes are necessary. In the absence of direct measurements of regurgitation, such as area and flow-rate measurements at the atrioventricular valve, if measurements of chamber volumes are not provided, the UKF has no information pertaining to the valve being regurgitant, and, consequently, a parameter estimation run yields a non-regurgitant valve where = EDV − ESV. In fact, even in the clinic regurgitant blood volume is assessed via the difference between EDV − ESV and . Consequently, it is critical to include measurements of ventricular volume along with measurements for (or alternatively ) to model regurgitation. While conductance catheters are the gold standard for measurement of time-varying ventricular volume information, this is rare in routine clinical practice. Commonly, the EDV and ESV are estimated from echocardiography, but more accurate measurements are possible from 3D echocardiography and MRI. These result in volumes at specific points of the cardiac cycle, such as end-diastole and end-systole, and not volume changes over the full cycle. Consequently, in addition to the methods employed for patient A, the EDV and ESV measurements are included in the UKF to provide information on valve regurgitation through the approach presented in §2.5. With these methods, the evolution of parameter estimate means for the 34 model parameters are shown in fig. 6 of [17]. The associated decrease in the variance of the parameter estimates is shown in figure 5. For this patient, all variances of the parameter estimates are reduced by at least a factor of four when the UKF is run for 40 cardiac cycles. Varying degrees of decrease in the parameter variances, however, are apparent, probably because of the difference in sensitivities. Owing to the inclusion of EDV and ESV measurements, the clinically measured regurgitation fraction of 25.3% is well captured in the model, which produces a regurgitation fraction of 22.7%. For further validation of the parameter estimates, through direct MRI measurements and indirect electrocardiogram and wedge pressure measurements, see [17] and the discussion in §3.2.2.

### 3.2. Inverse problems with heart rate variability

Results for heart rate variability are first presented for a synthetic dataset followed by a patient-specific case.

#### 3.2.1. Synthetic dataset

To demonstrate the efficacy of the proposed UKF method in taking measurements at different heart rates into account (see §2.7), a further reduced system of circulation is constructed. This system lumps all the systemic and pulmonary circulations into a single block and is shown in figure 6. For the heart activation functions, it is assumed that the ventricle activation duration and the ratio remain invariant with changes in heart rate. This implies that, at higher heart rates, the duration of ventricle activation remains constant while the decrease in the heart beat time period is reflected in a decrease in . Such an assumption is justified by clinical studies, for example see [48], which report that, at higher heart rates, the decrease in the systolic time interval is significantly smaller than the decrease in the diastolic time interval. Synthetic observations are generated by running the forward model with known parameter values at two different heart rates of 120 and 150 beats per minute (bpm). Pressure output is stored from the lower heart rate model and flow-rate output for the higher heart rate model. Small noise (2% of maximum values) is added to these two sets of outputs to generate measurements for parameter estimation. The combined set of measurements corresponds to []. The 16 parameters to be estimated are [ ].

Figure 7 shows how the two model states at different heart rates are corrected by the UKF. As expected, in each heart rate, the model state that is observed is corrected much faster than the unobserved state. For example, in the lower heart rate model, because only pressures are observed, the flow rates are corrected much more slowly than the pressures. Nonetheless, within 20 cardiac cycles (of the lower heart rate model), the UKF parameters converge very close to the true parameters. This is shown in figure 8. Finally, the evolution of UKF parameter variances for two cardiac cycles of the lower heart rate model are shown in figure 9. The results display a consistent decrease in the parameter variances with the observation of measurements. Here, within two cardiac cycles the reduction in the variance for all parameters is, at least, more than a factor of , with the exception of , where this factor is, at least, .

#### 3.2.2. Patient-specific case

In [17], it is noted that in patient B the pressure and flow-rate measurements were acquired at different heart rates of 106 and 140 bpm, respectively. In that study, the heart beat time period of the pressure measurements was artificially shortened to be consistent with the flow-rate measurements. Here, to account for this shortcoming, the method to deal with heart rate variability presented in §2.7 is applied to model the circulation of patient B. The aforementioned hypothesis that remains constant between different heart rates is made and is estimated as opposed to in [17]. The results of the UKF estimation are presented in figure 10. The evolution of estimate means and variances shows that all parameters are converged and that their variances decrease by a factor of at least 2^{5}. A comparison of model output (with UKF-estimated parameters) with the measurements is shown in figure 11. Comparing this figure with the results of [17] (reproduced in figure 11 in solid grey), it is observed that the measured venous flows, particularly and , are significantly better reproduced when heart rate variability is taken into account. It is also observed that the agreement between ventricular pressure in model output and the measurements is improved, particularly in relation to the duration of ventricular mechanical activation. The change of pressure–volume loop from lower heart rate (solid line) to higher heart rate (dashed line) is consistent with the assumption of and remaining constant when heart rate changes. It appears that at higher heart rates there is a restricted ventricular filling time, resulting in lower EDV. Consequently, there is lower ventricular preload and hence peak systolic pressure generation is reduced. Overall, due to better modelling of physics it is observed that consideration of heart rate variability improves the parameter estimates and the agreement between model output and the measurements, albeit at a higher computational cost (see §2.7). It should be noted that the force–frequency relationship [49], which suggests that cardiac muscles produce higher forces at higher frequencies, is not included in the model, and it is likely that inclusion of this phenomenon will result in further improvement of the results. This, however, is an area of future investigation.

Validation of the parameter estimation procedure can be performed in two ways: first, to consider clinical estimates of certain measurable parameters such as ventricular wall volumes (measured by MRI), valve annulus areas and mechanical activation times of heart chambers (qualitatively assessed by ECG readings); and second, by considering measurements of model output that were not employed in the parameter estimation procedure. The latter can include qualitative measures such as pulmonary venous wedge pressures and measures derived through echocardiography, and quantitative measures such as ventricular pressure and volume variations. For patient B, a twofold qualitative validation is performed: first, by the range of variation of pulmonary venous wedge pressure (measured to be versus produced by the model); and second by mechanical activation times of the single atrium and single ventricle which show an overlap between their contractions (implying that there is a significant time period when both the atrium and ventricle are contracting, which is indeed qualitatively observed in the ECG readings depicting the merging of the T and P waves). Quantitative measurement is also performed in two ways: first, by the measurement of ventricular pressure tracing that was not used for parameter estimation (see figure 11, row 2 from bottom to top, column 2, which shows this measurement in green and the model output in solid blue); and second by the measurements of end-systolic and end-diastolic ventricular volumes measured by MRI to be 13 ml and 30 ml, respectively, versus the model output of 12.7 ml and 30.3 ml, respectively.

When only a single heart rate is considered, see [17], with 20 parallel processes the UKF takes approximately 140 s for 200 discrete observations in one cardiac cycle. The same computation time when heart rate variability is considered is approximately 200 s.

While the proposed method works quite well, in certain cases it may not be possible to write a closed system of equations for the dynamical system without using a noisy measurement (for instance, as a forcing/boundary condition). This method may fail in such cases. The method may also not be applicable when only min/max/mean values are available as measurements without any full time-varying measurement at all. In such cases, it may be advantageous to consider other Bayesian methods for parameter estimation, such as MCMC; see for example [19]. Note that methods such as MCMC may also be used in general estimation problems, such as those studied in this article, albeit at a higher associated computational cost. Lastly, if the measurement noise is coloured (perhaps due to a correlated structure), then appropriate modifications, see for example [50], may be needed to the proposed method.

## 4. Conclusion

Several challenges of the UKF method for parameter estimation in haemodynamics problems are tackled. These include the prevention of Sigma point collapse to avoid ill-conditioning of the empirical covariance matrix, inclusion of min/max measurements in the method, and the choice of appropriate prior knowledge. Furthermore, an extension to the basic UKF method to assimilate measurements at different heart rates is presented. The methods are demonstrated on two patient-specific cases (one with and one without atrioventricular valve regurgitation). The efficacy of the proposed extension to the UKF method for assimilating measurements at different heart rates is first shown for a synthetic dataset as a proof of concept where the true parameters are known, and then on a patient-specific case with atrioventricular valve regurgitation. In the latter case, 34 parameters are estimated and the model predictions are shown to be significantly improved when heart rate variability is taken into account. The proposed methods enable analysis of patient-specific physiologies in a computationally efficient manner while taking measurement uncertainty into account. Lastly, these methods can be readily applied to other pathophysiologies.

## Authors' contributions

S.P. carried out the numerical experiments and drafted the manuscript. S.P., C.C., G.P. and I.E.V.-C. developed the mathematical models. S.P., G.P. and I.E.V.-C. conceived the study. C.B. and T.-Y.H. provided medical measurements, clinical input, and interpretation of results. All authors gave final approval for publication.

## Competing interests

We have no competing interests.

## Funding

This work was supported by the Leducq Foundation as part of the Transatlantic Network of Excellence for ‘Multi-scale modeling of single ventricle hearts for clinical decision support’ and a British Heart Foundation Clinical Research Fellowship FS/12/35/29566.

## Appendix A. Clinical data acquisition

All MRI acquisitions were performed on 1.5 T scanners (Philips Intera Achieva, Best, Netherlands; and Siemens Avanto, Siemens Medical Solutions, Erlangen, Germany). Functional, flow and 3D information was acquired in a routine pre-stage II clinical protocol. ECG-gated velocity-encoded phase contrast imaging sequences were used to acquire flow measurements in multiple locations under free breathing. Flow rates were calculated using the OsiriX open-source software (OsiriX Foundation, Geneva, Switzerland), through an in-house plugin. A contrast-enhanced 3D angiogram was acquired following administration of gadoteridol. Cardiac catheterization under general anaesthesia was performed following a routine clinical protocol in a bi-plane fluoroscopy suite (Siemens Medical Solutions USA Inc., Malvern, PA, USA). Pressure and haemodynamic measurements in several systemic and pulmonary arterial and venous locations were acquired through a fluid-filled catheter system. Valvular and ventricular function was assessed by routine echocardiography and pulse wave Doppler velocimetry.

- Received June 28, 2016.
- Accepted December 5, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.