Implementation and Calibration of a Deep Neural Network to Predict Parameters of Left Ventricular Systolic Function Based on Pulmonary and Systemic Arterial Pressure Signals

The evaluation of cardiac contractility by the assessment of the ventricular systolic elastance function is clinically challenging and cannot be easily obtained at the bedside. In this work, we present a framework characterizing left ventricular systolic function from clinically readily available data, including systemic and pulmonary arterial pressure signals. We implemented and calibrated a deep neural network (DNN) consisting of a multi-layer perceptron with 4 fully connected hidden layers and with 16 neurons per layer, which was trained with data obtained from a lumped model of the cardiovascular system modeling different levels of cardiac function. The lumped model included a function of circulatory autoregulation from carotid baroreceptors in pulsatile conditions. Inputs for the DNN were systemic and pulmonary arterial pressure curves. Outputs from the DNN were parameters of the lumped model characterizing left ventricular systolic function, especially end-systolic elastance. The DNN adequately performed and accurately recovered the relevant hemodynamic parameters with a mean relative error of less than 2%. Therefore, our framework can easily provide complex physiological parameters of cardiac contractility, which could lead to the development of invaluable tools for the clinical evaluation of patients with severe cardiac dysfunction.


INTRODUCTION
Heart failure corresponds to a clinical syndrome with a wide spectrum of symptoms, ranging from dyspnea and exercise intolerance to cardiogenic shock. It is caused by structural or functional cardiac abnormalities that result in low cardiac output, i.e., the inability of the heart to provide sufficient blood flow to satisfy the metabolic needs of the organism (Metra and Teerlink, 2017). It affects approximately 2% of the population, with a lifetime risk of developing heart failure of 20%, and a 5-years mortality of about 50% (Yancy et al., 2013). In the intensive care unit (ICU), cardiogenic shock represents 6% of admission, with an in-ICU mortality as high as 50% (Puymirat et al., 2017).
Evaluation of cardiac function is of crucial importance since it directly impacts treatment and therefore prognosis. In particular, the early and repeated evaluation of the left ventricular function is a key factor to guide therapies. This is particularly relevant in the ICU setting, where hemodynamic changes in unstable patients can occur rapidly (minutes). The aim of hemodynamic monitoring is therefore to accurately and rapidly identify changes in order to adapt treatment.
Hemodynamic assessment of cardiogenic shock is nowadays multimodal and includes clinical evaluation and paraclinical examination. Echocardiography for this situation is currently the cornerstone tool since it provides a complete evaluation of cardiac function and structure and can also predict response to treatment, in particular fluid-responsiveness (Vieillard- Baron et al., 2019). However, this examination is time-and resources-consuming (although new sensor technology may reduce the costs of these devices), hence evaluation based only on echocardiography may result in delays in management. Moreover, it cannot be used in a continuous manner. Other methods, either invasive or non-invasive, allow to indirectly assess cardiac output and volumes (Alhashemi et al., 2011;Monnet and Teboul, 2017;Nguyen and Squara, 2017).
Although these tools provide accurate estimates of ventricular function, evaluation of cardiac mechanics is fully characterized using the measure of pressure and volume of the ventricle during the entire cardiac cycle. From these values one can establish the ventricular pressure-volume diagram and determine the end-systolic and end-diastolic pressure-volume relationships. The former permit the determination of ventricular end-systolic elastance (Ees), a load-independent measurement of ventricular contractility (Burkhoff et al., 2005;Naeije and Manes, 2014). Nowadays these measurements require invasive techniques only available in the catheterization laboratory (Bonnet et al., 2017). Non-invasive methods have been developed, including magnetic resonance imaging (Bastos et al., 2019), or real-time three dimensional echocardiography (Seemann et al., 2019), but these cannot be easily implemented in the daily clinical practice.
To circumvent the drawbacks of these different methods, we aimed at developing a framework to characterize parameters of ventricular systolic function from easily accessible clinical data, namely systemic and pulmonary arterial pressures. The gold standard for such measurements is the use of invasively inserted intravascular catheters in the radial or femoral artery (systemic arterial pressure), or the pulmonary artery (pulmonary artery pressure). To obtain accurate pressure values, the measurements must be done at end-expiration and in the supine position, with standard zero reference at the level of the right atrium. Some variability of intravascular pressure measurement may still occur as a consequence of over-or underdamping of the signal (Romagnoli et al., 2014). We implemented and calibrated a deep neural network (DNN)-trained with data obtained from a lumped model of the cardiovascular system (Shi et al., 2011), as developed by Ursino (1998)-to model different levels of severity of left ventricular systolic dysfunction. This DNN takes, as input, systemic and pulmonary arterial pressure signals and returns, as output, the parameters of the lumped model that characterize left ventricular function, in particular Ees. The inputs are relevant measurable clinical values, whereas the predicted outputsas mentioned above-cannot be easily evaluated in clinical practice; these are in turn provided to the 0D model in order to recover all other clinical values, especially the ventricular pressure-volume curves.
The role of the proposed DNN is to solve an inverse problem which maps some of the outputs of the 0D model (i.e., systemic and pulmonary arterial pressures) to its underlying parameters. Of note, previous investigators, (Pennati et al., 1997;Pant et al., 2016;Schiavazzi et al., 2017) have proposed alternate methods for parameter identification of lumped parameter models in other physiological conditions, such as single-ventricle physiology with valve dysfunction (Pant et al., 2016). In this paper, we present a novel deep-learning-based approach to solve this issue.

Global Framework
In this section, we focus on the online pipeline of our method (i.e., the prediction procedure) which is depicted in Figure 1, whereas the training process is described in section 2.3 for general DNNs and section 2.4 for our specific application. Input parameters were systemic and pulmonary arterial pressure curves. Theses curves are represented in the frequency domain by their Fourier coefficients. As explained in more detail in section 2.4, the Fourier coefficients were subsequently arranged in a three-dimensional tensor which served as input for the DNN. The output of the DNN were the four parameters characterizing left ventricular function in the 0D model, namely E max,lv , E max lv,0 , G E max,lv , and k E,lv (see section 2.2 for more details). In the last step, the predicted parameters were provided to the 0D model to simulate and recover all the other values.

Lumped Parameters Model of the Cardiovascular System
The lumped model of the cardiovascular system has been extensively described in the original paper by Ursino (1998). The model provides a mathematical description of the entire cardiovascular system with time-varying elastance of the left and right ventricles. It also includes the afferent carotid baroreceptor pathway, the sympathetic and vagal efferent activities, the splanchnic and extrasplanchnic systemic circulation, and the pulmonary circulation. A representation of the model is presented at the bottom of Figure 1. For each of the vascular icompartment, mass and momentum conservation are given by, respectively: and where V i is the volume of the i-compartment, Q in and Q out the inflow and outflow rate, P in and P out the inlet and outlet pressure, R i the resistance, and L i the inertance. The pressure and flow relationship reads Inertance is taken into account for arterial compartments where acceleration of blood is significant, while neglected for the other ones. Valve opening is determined by its pressure gradient across the valve, i.e., it opens when P in > P out . Left and right ventricle contractility is modeled as a time-varying elastance. Efferent pathways act on several parameters, namely heart period, left and right ventricle contractility, unstressed volumes (defined as the volume when pressure is equal to zero), and both splanchnic and extrasplanchnic peripheral artery resistance. Parameters affected by autoregulation are surrounded by a gray line in Figure 1.
In a previous work (Bonnemain et al., 2013), we modified the model to account for left ventricular systolic failure at different levels of severity and validated it with clinical data. To reproduce various degrees of left ventricular systolic failure, we varied the following parameters: • E max,lv [mmHg/ml], reference value of the ventricle elastance, • E max lv,0 [mmHg/ml], reference value of the ventricle elastance in absence of autoregulation, • G E max,lv [mmHg/ml/(spikes/s)], maximum baroreceptor gain, • k E,lv [1/ml], steepness of the pressure-volume curve.

Figure 2
shows left ventricular pressure-volume loops at different degrees of severity of heart failure. All other haemodynamic values are presented in Bonnemain et al. (2013). For each compartment, pressure, volume and flow can be extracted; it is therefore straight-forward to retrieve the ventricular pressurevolume curves.

Deep Neural Network
The term machine learning is typically used to indicate a family of algorithms and methods which, broadly speaking, are capable of identifying structures and patterns in complex and (usually) high-dimensional information. Due to the abundance of data obtained during the last decade and the ever increasing computational power of processors, graphics processing units (GPUs) and supercomputers, these methods have gained particular attention in all the fields of medicine (Komorowski et al., 2018;McKinney et al., 2020).
DNNs are a class of data-hungry machine learning algorithms notably suited for classification and regression tasks. They are used to approximate unknown mappings f(·) of the form y = f(x; ω), where x and y are input and output of the model, respectively, and ω indicates a set of parameters of the network. The function f(·), which, in reality, can be as general as possible, is represented in DNNs as a composition of easily computable parameterized functions f 1 , f 2 , . . . , f l (where l is the number of layers of the network), i.e., f = f l • . . . • f 2 • f 1 . In this work, we focused on the most basic type of DNNs, i.e., multilayer perceptrons, which, despite their simple structure, have been mathematically proven to satisfy desirable approximation properties. In particular, multi-layer perceptrons with one hidden layer or more are universal approximators, see e.g., (Cybenko, 1988;Mhaskar, 1993). In contrast to more advanced architectures (e.g., convolutional neural networks), multi-layer perceptrons are solely based on a sequence of fully-connected layers of the form where the subscript i refers to the ith layer of the network, x i = y i−1 is the corresponding input, W i and b i are a weight matrix and a bias vector, and σ i is a non-linear activation function. Among the most common activation functions are e.g., the rectifier linear unit ReLU(x) = max(0, x), sigmoid(x) = 1/(1 + exp (−x)), and tanh(x). In this paper, we consider the ReLU and sigmoid activation functions for the hidden (i.e., f 1 , . . . , f l−1 ) and output (i.e., f l ) layers, respectively. The weight matrix W i and the bias vector b i represent the parameters ω i of the ith layer. In order to find an appropriate numerical value of the parameters ω i , for each i = 1, 2, . . . , l, the DNN has to be trained on a set of input-output pairs Q = {{x 1 , y 1 }, {x 2 , y 2 }, . . . , {x n t , y n t }} which constitute the so called training dataset. The ultimate goal of the training process is to find ω = {ω 1 , ω 2 , . . . , ω l } such that the DNN f(·, ω) performs a "good" fitting on Q, i.e., f(x k , ω) ≈ y k for every k = 1, 2, . . . , n t . This is typically done by introducing a loss function acting on the training dataset L ω (Q) and by minimizing it by means of optimization algorithms (usually, gradient or stochastic gradient descent). For a complete review of DNNs and their optimization, we refer the reader to Goodfellow et al. (2016).

The Offline Phase: Data Generation and Training of the Deep Neural Network
In the data generation phase, the 0D model described in section 2.2 is solved N s times by randomly sampling E max,lv , E max lv,0 , G E max,lv , and k E,lv . Specifically, we determined for each of these parameters a range of physiological or pathological (severe heart failure) values; for details about these ranges see (Bonnemain et al., 2013). Upper and lower bounds are reported in Figure 2. Each training sample was generated from the 0D model, by assigning a random individual value (drawn from a uniform distribution defined on the corresponding admissible range) for each of the four parameters of interest, and by keeping constant all the other parameters of the 0D model. Other parameters include values of resistance, compliance, unstressed volume, inertance, as well as values describing heart function, and values related to autoregulation, as presented in detail by Ursino (1998). This approach allows generating a dataset which includes a rich representation of all possible pathological conditions leading to left heart failure. Regarding sampling strategies, owing to the negligible computational cost required the generate each random sample, we did not consider alternative sampling strategies (e.g., Latin Hypercube or orthogonal sampling). The output of the 0D model comprises multiple time dependent one-dimensional signals. We focused only on quantities that are easy to measure in the clinical practice, such as systemic arterial pressure and pulmonary arterial pressure. Let us denote (0, T) the time interval over which the simulation of the 0D model is performed, t 1 = 0, . . . , t n t = T a set of timesteps such that the timestep size t = t m+1 − t m is constant for all m = 1, . . . , n t , and u k,i (t m ) the value of the ith quantity of interest (e.g., systemic or pulmonary arterial pressure) of the kth sample at timestep t m . In order to better capture the time component of the quantity of interest, we decided to operate on its representation in the frequency domain. More precisely, we computed the Discrete Fast Fourier Transform (DFFT) of the signal u k,i (t m ) and we recovered the coefficients a k,i j and b k,i j such that, for all k = 0, . . . , n t , where M = n t /2 and ω = 2π/T. We note that Equation (5) must be slightly modified to account for the case of n t odd; we refer to Quarteroni et al. (2014) for details. We introduce for sake of clarity the vector c k,i = [a k,i 0 , . . . , a k,i M , b k,i 1 , . . . , b k,i M ], i.e., the vector containing all the 2M + 1 Fourier coefficients. In this work we selected the first 5% of Fourier coefficients, i.e., M = 5. Cutting higher modes allows to eliminate noise without loosing information. We studied the influence of the number of Fourier coefficients used for the DNN input on the DNN performance. Data from these studies are provided in Supplementary Data Sheet 2. Overall, results indicate that the best accuracy of the DNN is obtained using M = 5 (11 Fourier coefficients).
The training dataset Q for the DNN is composed of pairs {x k , y k }, k = 1, . . . , N s , where x k = [c k,1 , . . . , c k,N i ], N i being the number of desired input signals selected from the output of the 0D model, and y k = [E k max,lv , E k max lv,0 , G k E max,lv , k k E,lv ]. In our implementation, the input data is organized in a three-dimensional tensor X lmn = c l,m n with dimensions N s × N i × (2M + 1) as depicted in Figure 1. Therefore, the forward pass of the DNN consists in mapping the frequency coefficients of a handful of time dependent signalsin our case, systemic and pulmonary arterial pressure-to the corresponding values of the hidden parameters of the 0D model (as already mentioned in section 1, the DNN aims at solving an inverse problem). Regarding the choice of loss function L ω (Q), in this work we focused on the mean squared error MSE ω (Q) = 1/N s N s k=1 |y k − f(x k ; ω)| 2 , which is well-suited to regression tasks.

RESULTS AND DISCUSSION
The 0D model was implemented in the Modelica language 1 , a non-proprietary, object-oriented, equation based language which conveniently models complex physical systems. We used OpenModelica 2 as a modeling and simulation environment. The DNN was implemented using TensorFlow 3 . Discrete fast Fourier transform computation and statistical analysis were performed with Matlab 4 .
As described in section 2.4, N s = 10 ′ 000 samples were generated from the 0D model; 5% of these were reserved to the test set, which was employed to evaluate the performance of the trained network. Of the 95% samples remaining, 80% represented the train set and 20% the validation set. After considering several DNN architectures, we opted for a multi-layer perceptron with 4 fully connected hidden layers featuring 16 neurons per layer. Our choice was driven by the value of loss function achieved during training and by the absence of overfitting (which occurs whenever the error on the train set is considerably smaller than that on the validation set); in cases of architectures achieving similar performance, we chose the network composed of the smallest number of neurons. All the data regarding the considered DNN architectures and the relative performances can be found in the Supplementary Data Sheet 1. The learning rate is set at a value of 0.001. The Adam optimization algorithm was used to update the network weights during training. In Figure 3 we report the DNN performance for the 4 predicted parameters on the test set. Each graph represents the plotting of values of real parameter (x-axis) and the predicted parameter (y-axis). We note that the optimal prediction corresponds to the bisector. The architecture shows good accuracy especially for E max,lv , E max lv,0 , and k E,lv . In contrast, the prediction of G E max,lv appeared less accurate, which suggests that the 0D model is significantly less sensitive to this parameter.
To support this hypothesis, we performed a sensitivity analysis of the output of the 0D model with respect to small variations of the 4 parameters. The mean value of parameter P is given by P mean = (P min + P max )/2, where P min and P max are the lower and upper bounds depicted in Figure 2. For each parameter P, we run a set of 100 simulations by assigning the value P = P mean (1 + ǫN(0, 1)), where ǫ = 0.01 and N(0, 1) is a normal distribution with 0 mean and a standard deviation equal to 1, while keeping the other 3 parameters constant and equal to their respective mean values. The results are reported in Table 2. The output of interest shows a noticeably lower standard deviation with respect to small variations of G E max,lv , highlighting a lower sensibility of the 0D model for this parameter.
As depicted in Figure 1 and described in section 2.1, the main objective of our study was to recover a wealth of haemodynamic values from systemic and pulmonary arterial pressure signals. Table 1 shows the results of the 0D model simulations solved for the real and predicted parameters. Specifically, the 0D model was run for each 4-parameters set of the test set (5% of total samples, i.e., 500 samples), both with the real and the predicted parameters. Relevant haemodynamic values were extracted and compared. These values were: systemic and pulmonary pressures (mean, systolic, and diastolic), heart rate, left ventricular ejection fraction, left ventricular end-diastolic and end-systolic volumes, cardiac index, and pulmonary capillary pulmonary wedge pressure. Minimal, maximal, mean, and standard deviation of the real and predicted values were calculated. The calculated difference between real and predicted values was defined as the error, for which standard deviation with 95% confidence interval were computed. All values are reported in Table 1.
Our test set included haemodynamic values spanning a range of clinical situation from normal physiology to severe heart failure. Overall our results indicate an excellent correlation between real and predicted values, with mean relative error never superior to 2%, and a narrow 95%-confidence interval. The accuracy of this model is relevant to the clinical setting, where these hemodynamic variables are susceptible to show wide variations.
Furthermore, the small error in our model should be regarded as perfectly acceptable, when considering the variability of intravascular pressure measurements in the clinical setting, notably related to underdamping and resonance phenomena (Romagnoli et al., 2014). The same is true for clinical measurements of ventricular volumes (Bastos et al., 2019). With respect to ventricular ejection fraction, it is worth mentioning that its measurement by standard modern techniques (cardiac computerized tomography, radionucleotide and invasive ventriculography, echocardiography or magnetic resonance imaging) presents a bias of less than 5% (Pickett et al., 2015), which reinforces the relevance of our model, which leads to a relative error smaller than 2%. Finally, new devices and noninvasive approaches for continuous pressure measurement may broaden the applicability of the presented framework (Proena et al., 2016), while new methods for haemodynamic monitoring could provide innovative perspectives (Pour-Ghaz et al., 2019).
Our model appears slightly better to predict pressures than other parameters, which may reflect the fact that the DNN takes pressure, but not volumes values, as input. In addition, with respect of heart rate, it can be noticed that minimal predicted and real values were identical (60/min), which simply reflects the lowest boundary of this parameter in the 0D model.

LIMITATIONS
We acknowledge the three following limitations in our study. First, our model is adapted to the pathophysiology of left heart failure, and therefore cannot apply to other condition of cardiac failure, such as right heart dysfunction. Secondly, we choose to focus only on 4 parameters of interest. Therefore, it is likely that including additional parameters (for instance a resistance parameter associated with vascular dysfunction) could have affected the resulting accuracy of the network. Third, even though the 0D model includes a component of autoregulation, it does not take into account respiratory physiological variables, which may significantly influence cardiac function (concept of heart-lung interactions), notably in the context of mechanical ventilation, which is frequently provided in patients with severe cardiac failure. Fourth, our study used exclusively synthetic data. Using real patients data (possibly affected by noise) to validate our framework will be essential for its applicability in the clinical setting, and such validation will be the matter of future investigations. Furthermore, we cannot rule out that using real clinical data could affect the ability to train the proposed network, an issue that will require additional studies.

CONCLUSIONS AND PERSPECTIVES
In this work we presented a framework that allows to predict, from physiological data (systemic and pulmonary arterial pressure signals), various parameters of left ventricular function and other haemodynamic variables. The 0D model used to train the DNN in our study allows to describe a wide spectrum of pathophysiological alterations pertaining to left heart failure. In turn, the DNN displayed remarkable accuracy to recover relevant haemodynamic parameters. Our study represents a first step toward the development of automated tools providing helpful information on heart function. For instance, the application of our framework could assist in the real time detection of left ventricular dysfunction, especially in ICU patients, who are subject to continuous monitoring of intravascular pressures. Also, our tool could be useful to assess the left ventricular function during mechanical circulatory support, such as left ventricular assist device (LVAD) or veno-arterial extracorporeal membrane oxygenation (ECMO).
Future studies, that will validate the model with patient related data, will be necessary for the clinical implementation of our tool. Moreover, our framework will be of great interest for the future extensions of the 0D model to other circulatory pathologies or heart disease as well as to mechanical circulatory support.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation. Implementation of the lumped parameter model can be found at: https://bitbucket.org/jbonnema/lumped_parameter_ model/src/master/.

AUTHOR CONTRIBUTIONS
JB contributed to the design and the implementation of the research, to the analysis of results and to the writing of the manuscript. LP contributed to the design and the implementation of the research, and to the writing of the manuscript. SD contributed to the design of the research, to the analysis of the results, and to the writing of the manuscript. LL contributed to the analysis of the research and to the writing of the manuscript. All authors contributed to the article and approved the submitted version.