Identification of cardiac afterload dynamics from data

The prospect of ex vivo functional evaluation of donor hearts is considered. Particularly, the dynamics of a synthetic cardiac afterload model are compared to those of normal physiology. A method for identification of continuous-time transfer functions from sampled data is developed and verified against results from the literature. The method relies on exact gradients and Hessians obtained through automatic differentiation. This also enables straightforward sensitivity analyses. Such analyses reveal that the 4-element Windkessel model is not practically identifiable from representative data while the 3-element model underfits the data. Pressure–volume (PV) loops are therefore suggested as an alternative for comparing afterload dynamics.


Ex vivo heart evaluation
Today, donor hearts are routinely discarded due to uncertainty about their function. A method of functional evaluation of donor hearts therefore has the potential to increase the availability of heart transplantation.
Such ex vivo (outside of the body) functional evaluation requires a system that mimics vascular dynamics-the afterload-and can be controlled to emulate a broad range of recipient physiology and working conditions. The afterload can be decoupled from heart dynamics by simultaneously measuring aortic flow and pressure, and relating them through an impedance model. Here, we delimit our focus to Windkessel models, being a class of linear time-invariant (LTI) arterial impedance models introduced by Otto Frank in the late 19 th century and commonly employed both academically and clinically.
We propose a method for identifying the Windkessel models from data, and perform sensitivity analyses to quantify uncertainty of the identified models. Three data sources are considered: previously published in vivo human data; porcine heart beating in vivo; porcine heart beating ex vivo against a synthetic afterload. The data and code used to generate the results are available on GitHub, see Pigot (2021).
This work focuses on technical aspects related to identifiability and parameter identification; details on the device and data generating experiments will be presented elsewhere.

The arterial Windkessel
Several variants of the Windkessel model have been proposed in the literature, see for instance Westerhof et al. (2009). Here we have employed the parallel 4-element Windkessel model, being one of the more general formulations. It can be expressed in terms of the circuit analogy shown in Fig. 1, that relates pressure (potential) to antegarde flow (current) through a dynamic impedance , where we use subscript to denote continuous time. The passive component parameters are described in table 1.
Using the circuit analogy of Fig. 1 we next derive the expression for the transfer function ( ). In the Laplace domain the pressure (voltage) across resistance , inertance (inductance) and compliance (capacitance) relate to the flow (current) flowing though each element through = , = , and = , respectively. Denoting by the (pressure) potential between the resistances according to Fig. 1, Kirchhoff's current law yields From (1) we can eliminate = /(1 + ) to obtain the transfer function The Windkessel model (2) is hence a parallel interconnection between and two first-order systems. Note that the characteristic resistance only impedes accelerating flows, since for steady flows, the impedance of the inertance is 0. This explains why the static gain of (2) is (0) = (and not, as some might assume + ). The 4-element Windkessel also comes in a less widespread series form, used in for example Gellner et al. (2020). In the series form, the inductance is connected in series with , and its transfer function is obtained by replacing the last term in (2) by . The static gain is + , but the transfer function is improper, as the term corresponds to an unfiltered differentiation of the input.
Introducing the state = [ 1 2 ] , defined through results in the state space form S 4 with matrices defined through with 1 being the volume (charge) in . From (4) it is directly visible that the system has two poles corresponding to time constants 1 = and 2 = / .

Experiments
The porcine data used in this work were recorded from two 65 kg Swedish pigs (sus scrofa domesticus). Large-animal experiments were needed to obtain a reliable characterization of the synthetic afterload module under consideration. All institutional and national guidelines for the care and use of laboratory animals were followed and approved by the appropriate institutional committees. The animals were treated in compliance with Directive 2010/63/EU, The European Parliament (2010). The study ran under ethics permission M174-15, issued by "Malmö/Lunds Djurförsöksetiska Nämnd" (local REB).
In both the in vivo and ex vivo experiments, a CardioMed CM4000 ultrasonic transit time flow meter (Medistim ASA, Oslo, Norway) was secured to the ascending aorta. Aortic pressure was measured using a Meritrans DTXPlus pressure transducer (Merit Medical, Singapore). In-house developed data acquisition hardware and software were used to log the signals at 200 Hz.
Different hearts were used for the in vivo and ex vivo data. The in vivo data reflects normal sinus rhythm in rest. The ex vivo data was recorded after 24 hours of cold non-ischemic perfusion, as described in Steen et al. (2016), beating against an actively controlled synthetic afterload akin to the device described by Gellner et al. (2020).

Model formulation
We relate the continuous-time flow ( ) to the aortic pressure ( ) through = * + , where is the impulse response of a an LTI system model with transfer function , and is the output error signal of the model, also referred to as the residual.
We do not have access to and directly, but only to corresponding measurement time series , , each of elements. The time series will here be considered equitemporarily sampled at a period of ℎ, although the proposed methodology is readily applicable also under irregular sampling schemes.
Applying the zero-order-hold operator X ℎ we thus obtain [ , , , ] = X ℎ [ , , , ] that relate through = * + and discrete time state space realization S : In order to evaluate the residual, , we simulate S with as input.The residual can be decomposed into one modelmismatch term and one noise term. The former typically arises from a candidate under-modelling the data , . The latter arises if cannot be fully explained by . This is for instance the case if and have been subjected to measurement noise. In this work we do not discern between the two contributors to .
To determine the initial state vector for this simulation, we could choose an arbitrary value, e.g. 0 = 0 and drive S with a repeated stack [ . . . ] with sufficiently many repetitions of for the transient caused by 0 to fade, and then discard all but the last simulated output samples. An efficient and approximation-free alternative is to directly enforce the corresponding "periodic stationarity" condition 0 = . Simulating the system forward in time over one cardiac cycle then gives Using that = 0 we solve (5) for where, for numeric robustness, the sum in (5) is preferably computed as a convolution, rather than term-by-term.

Parameter identification
As stated in Sec. 2.2 we do not explicitly consider observation model (or other noise source) stochastics in this work, and therefore proceed with output error identification. The objective is hence to minimize the residual cost to identify = arg min 0 ( ).
The solution of (8) will also minimize = 1 2ℎ 2 2 , under the assumption that is piece-wise constant between the samples in , and that the signals and have been subjected to analog filtering effectively removing spectral power above the Nyquist frequency (2ℎ) −1 . The latter holds true for the data considered here. The former constitutes a valid approximation as long as 1/ℎ is large compared to the magnitudes of the poles of . Validity of this approximation has been ensured through sampling at a high frequency compared to the cardiac cycle dynamics.
The optimization problem (8) is generally non-convex in . We therefore employ a multiple initialization procedure. Being the static gain of the system, is initialized with the mean output (pressure) divided by mean input (flow). The initial values of the remaining parameters are drawn from a multivariate uniform distribution, and we retrospectively verify that this distribution covers the identified parameters .
For each initialization point, Newton's method is then used to obtain a cost minimizer candidate, and the one of them that minimizes gets denoted .
A dual number implementation (see Revels et al. (2016)) of the zero-order-hold operator X ℎ enabled exact (down to machine precision) forward-mode automatic differentiation and thus evaluation of the gradient ∇ ( ) and Hessian ∇ 2 ( ) required in each Newton iteration. This enables us to identify the continuous-time model parameter directly, without the need of finite difference, or other, approximations. It also means that we can obtain an exact evaluation of the Hessian = ∇ 2 ( ) of the cost with respect to the parameter , evaluated at the optimum .

Persistence of excitation
A classical result for identification of LTI systems is that the order of persistence of excitation (PE) of an input sequence determines whether is sufficiently informative to distinguish parameter candidates and given some model class, or structure, . Particularly, PE of order at least is required for to be sufficiently informative to distinguish any pair of transfer functions of parameters. See for example Mareels et al. (1987) for further detail.
The PE degree can be determined in several (equivalent) ways, one being through a rank condition on the expectation of the auto-correlation matrix of : where ( ) is expectation of the auto-correlation of , with respect to some stochastic observation model relating to . In absence of stationary (additive) noise model, it is customary to assume that the observation ( ℎ) is an unbiased and consistent estimate of ( ) and use the observed auto-correlation, defined through (9) with where˜ ( ) is typically defined as ( ) for 1 ≤ ≤ , and 0 otherwise. Since we are dealing with signals of periodic nature here, we will instead use the periodic expansioñ ( ) = (( mod ) + 1). We can efficiently evaluate the corresponding circular sample auto-correlation sequence where is the × DFT matrix, its conjugate, and denotes element-wise matrix-matrix multiplication.
The rank condition states that in order for to be PE of order at least , the corresponding Φ needs to be positive definite. A major concern here is that any can be turned into a signal of arbitrarily high PE order by adding a signal of independent samples ∼ N (0, 2 ), even with arbitrarily small > 0. From a practical perspective it is therefore more relevant to consider the spectrum of Φ .

Sensitvity analysis
Next, we investigate how much model fit, expressed in terms of the cost , deteriorates if the optimal parameter is subject to an additive perturbation with magnitude 2 . The Taylor series expansion of the perturbed cost is where the residual ( ) is a linear combination of monomials in the components of , each with degree at least 3. If 2 is small, then the contribution of ( ) to is small. In that case the cost increment is well approximated by where the Hessian = ∇ 2 ( ) is a symmetric real matrix and therefore has a singular value decomposition according to (13). The singular vectors make up the columns of the unitary matrix = [ 1 .
. . ] and we can write the quadratic form as For a fixed 2 , (14) is minimized (maximized) when is parallel to the singular vector corresponding to the smallest (largest) singular value . This reveals in what direction a small move away from contributes least (most) to increase . Further, the fraction between the largest and smallest singular value, being the condition number of , reveals the relative change in cost when moving a small (infinitesimal) distance 2 in the least and most sensitive directions, respectively. A large condition number is therefore an indicator of over-parametrization with respect to the experimental data.

Identified models
The identification method was first validated against a previously published human data set-henceforth refered to as the human data-by digitizing the waveforms in Figure 4 of Stergiopulos et al. (1999) using Webplotdigitizer, Rohatgi (2020). Numerical values are reported in table 2; time domain model outputs and corresponding pressurevolume (PV) loops in Fig. 2; output error residuals in Fig. 5. As seen in Fig. 2, the parameter values identified using the method of Sec. 2 reproduce the appearance of the results in Stergiopulos et al. (1999). For reasons explained above in Sec. 2.5, 2-and 3-element Windkessel models were also identified. Bode plots of the identified models are shown in Fig. 6. These indicate that for the experimental data an almost perfect pole-zero cancellation takes place in the 4-element models, as can be verified by inserting the parameter values of table 2 into (2). This also explains why corresponding time-domain plots are visually indistinguishable.

Persistence of excitation
The 10 largest singular values of the auto-correlation matrices Φ for the inputs of the human, in vivo, and ex vivo data sets are shown in Fig. 7. According to Sec. 2.4, the figure indicates that identification of 2 to 4 parameters may be feasible, which prompts further consideration of parameter sensitivity in the fitted models. Since the input is generated by the heart, it cannot be arbitrarily changed to increase excitation and elucidate more parameters.
One can note that two parameters (degrees of freedom) are necessary to reproduce arbitrary diastolic and systolic  pressure levels. For the 4-element Windkessel model S 4 , this leaves two parameters to influence the shape of the model output; one "shape" parameter for S 3 ; none for S 2 .

Parameter sensitivities
The singular value decompositions = Σ of the Hessian = ∇ 2 ( ), evaluated at the identified parameter = are shown in table 3. The last column of the matrices indicate that the least certain direction in parameter space coincides with for the human data. The in vivo and ex vivo data instead result in uncertain estimates of . Recalling from Sec. 1.2 how limit cases of and correspond to the 2-and 3-element Windkessel structures, it is not surprising-given the indicated sensitivities-that the 3-element models explain the experimental data almost identically well, while the 2-element counterparts fail to do so based on lacking degrees of freedom, as mentioned in Sec. 3.2.

DISCUSSION
A method for identification of continuous time dynamics from time series data has been proposed with the objective of comparing the dynamics of a synthetic afterload with those of normal physiology. Upon validation against previously published results from Stergiopulos et al. (1999), the method was therefore applied to two porcine data sets: one collected in vivo, featuring normal physiology; one collected ex vivo, with the heart working against a synthetic afterload, as may be used in functional evaluation of donor hearts prior to possible implantation.
We can note that model fit of the human data is slightly better than for our experimental data. This is presumably caused by how our measurement setup was devised, and we are planing to conduct similar but refined measurements to investigate this. Nonetheless, the output error residuals from all three data sets, shown in Fig. 5, suggest that the Windkessel structure under-models representative pressure-volume data. At the same time, the persistence of

CONCLUSION
The well-established 4-element parallel Windkessel model is not reliably identifiable from representative pressurevolume time series data. Particularly, the peripheral resistance , being the static gain of the model, is reliably identifiable, while inertance parameter is unidentifiable.
For the sake of comparing afterload impedance dynamics, it is therefore advisable to look directly at the PV loop representation, rather than comparing Windkessel model parameters.

ACKNOWLEDGMENTS
This work was partially funded by the Swedish Research Council (grant 2017-04989) and the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation. The authors from the Department of Automatic Control are members of the ELLIIT Strategic Research Area at Lund University.