Rapid Estimation of Left Ventricular Contractility with a Physics-Informed Neural Network Inverse Modeling Approach

Abstract

Physics-based computer models based on numerical solution of the governing equations generally cannot make rapid predictions, which in turn, limits their applications in the clinic. To address this issue, we developed a physics-informed neural network (PINN) model that encodes the physics of a closed-loop blood circulation system embedding a left ventricle (LV). The PINN model is trained to satisfy a system of ordinary differential equations (ODEs) associated with a lumped parameter description of the circulatory system. The model predictions have a maximum error of less than $5\%$ when compared to those obtained by solving the ODEs numerically. An inverse modeling approach using the PINN model is also developed to rapidly estimate model parameters (in $\sim$ 3 mins) from single-beat LV pressure and volume waveforms. Using synthetic LV pressure and volume waveforms generated by the PINN model with different model parameter values, we show that the inverse modeling approach can recover the corresponding ground truth values, which suggests that the model parameters are unique. The PINN inverse modeling approach is then applied to estimate LV contractility indexed by the end-systolic elastance $E_{es}$ using waveforms acquired from 11 swine models, including waveforms acquired before and after administration of dobutamine (an inotropic agent) in 3 animals. The estimated $E_{es}$ is about 58% to 284% higher for the data associated with dobutamine compared to those without, which implies that this approach can be used to estimate LV contractility using single-beat measurements. The PINN inverse modeling can potentially be used in the clinic to simultaneously estimate LV contractility and other physiological parameters from single-beat measurements.

Keywords: Cardiac contractility, Patient-specific Modeling, Physics-informed Neural Network, Lumped Parameter Model, Parameter Estimation, Sensitivity Analysis

1 Introduction

Computational and/or mathematical models are increasingly used as patient-specific tools to understand cardiovascular diseases and treatments [1, 2, 3, 4] as well as to quantify changes in heart function (e.g., cardiac contractility) [5]. These models are usually developed based on numerical methods (e.g., finite element method) that are used to solve the partial differential equations (PDEs) and/or ordinary differential equations (ODEs) governing the key driving physics. Such models are often referred to as “physics-based” models [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]. Personalizing a physics-based model for each patient, however, requires estimating its parameters through inverse modeling, a process that often necessitates multiple runs of the numerical model. This process is time-consuming and highly computationally expensive, making it impractical when dealing with a large number of patient datasets or in time-sensitive clinical situations [17, 18].

Machine learning (ML) models have emerged as a promising solution to issues related to the lengthy computational time associated with patient-specific modeling [19] and modeling of biological systems [20, 21, 22, 23, 24, 25, 26]. While most ML models are trained purely with substantial amount of labeled data (i.e., data-driven ML model), another class of ML models that are directly trained to satisfy the governing equations of the operating physics (i.e., physics-based ML model) is increasingly being developed. Specifically, physics-informed neural networks (PINNs) where residuals of the governing equations are incorporated in the loss function during training with a small set of data or without data are increasingly developed [27] and applied in various biomechanical studies [28, 29, 30, 31, 32, 33, 34, 18, 35]. To the best of our knowledge, a PINN model describing hemodynamics in the cardiovascular system with varying cardiac function has not been developed.

Here, we developed a PINN modeling framework for rapid prediction of hemodynamics in the closed-loop cardiovascular system that is described by a system of ODEs. We use the PINN model to perform global sensitivity analysis that would otherwise be computationally prohibitive if numerical methods were used to solve the ODEs. We also apply the PINN framework in an inverse modeling approach to rapidly estimate physiologically-relevant parameters such as the left ventricle contractility from animal-specific hemodynamics waveforms. We demonstrate the translational potential of the PINN inverse modeling approach as a clinical tool for quantifying heart function by showing that this approach can capture an increase in left ventricle contractility with hemodynamics waveforms acquired after the administration of dobutamine.

2 Methods

2.1 Mathematical model

2.1.1 Blood circulation

The closed-loop lumped parameter circulatory model is comprised of five compartments (Fig. 1). The volume change of each compartment is governed by inflow and outflow rates as described by the following system of ODEs

Figure 1: The electrical equivalent diagram of the closed-loop blood circulation. lv, left ventricle; ao, aorta; art, peripheral artery; vc, vena cava; la, left atrium; av, aortic valve; mv, mitral valve.

$$\frac{dV_{lv}}{dt} =q_{mv}-q_{av}\,,$$ (1a)

$$\frac{dV_{ao}}{dt} =q_{av}-q_{ao}\,,$$ (1b)

$$\frac{dV_{art}}{dt} =q_{ao}-q_{art}\,,$$ (1c)

$$\frac{dV_{vc}}{dt} =q_{art}-q_{vc}\,,$$ (1d)

$$\frac{dV_{la}}{dt} =q_{vc}-q_{mv}\,,$$ (1e)

where $V_{lv}$, $V_{ao}$, $V_{art}$, $V_{vc}$, $V_{la}$ are volumes of the left ventricle (lv), aorta (ao), peripheral arteries (art), vena cava (vc), and left atrium (la), respectively. Flow rates in the aorta $q_{ao}$, peripheral arteries $q_{art}$ and vena cava $q_{vc}$, as well as across the mitral valve $q_{mv}$ and aortic valve $q_{av}$ are determined from their corresponding resistances (i.e., $R_{ao}$, $R_{art}$, $R_{vc}$, $R_{mv}$, $R_{av}$) and pressure differences between the connected storage compartments by

$$q_{ao} = \dfrac{P_{ao}-P_{art}}{R_{ao}}\,,$$ (2a)

$$q_{art} = \dfrac{P_{art}-P_{vc}}{R_{art}}\,,$$ (2b)

$$q_{vc} = \dfrac{P_{vc}-P_{la}}{R_{vc}}\,,$$ (2c)

$$q_{av} = \dfrac{\max(P_{lv}-P_{ao},\,0)}{R_{av}}\,,$$ (2d)

$$q_{mv} = \dfrac{\max(P_{la}-P_{lv},\,0)}{R_{mv}}\,.$$ (2e)

Note that Eqs. (2d) and (2e) are piecewise equations describing the flow rate by the valves’ opening and closing. These equations do not exhibit smoothness in their derivatives, which can pose an issue when training the PINN model. To circumvent this issue, a smooth approximation of the maximum function given by

$$\max(x-y,\,0)\approx\dfrac{1}{\alpha}\text{softplus}\left(\alpha(x-y)\right)\,,$$ (3)

is used instead, where

$$\text{softplus}(\alpha(x-y))=\begin{cases}\log(1+\exp(\alpha(x-y)));&\alpha(x-y)\leq 20\\ \alpha(x-y);&\alpha(x-y)>20\,,\end{cases}$$ (4)

and $\alpha$ is a parameter that controls the smoothness of the function approximation. Larger values of $\alpha$ lead to a closer approximation to the actual maximum function. Pressure in the blood vessels is given by

$$P_{ao} = \dfrac{V_{ao}-V_{ao,\,r}}{C_{ao}}\,,$$ (5a)

$$P_{art} = \dfrac{V_{art}-V_{art,\,r}}{C_{art}}\,,$$ (5b)

$$P_{vc} = \dfrac{V_{vc}-V_{vc,\,r}}{C_{vc}}\,,$$ (5c)

where $V_{ao,\,r}$, $V_{art,\,r}$, and $V_{vc,\,r}$ are the prescribed resting volumes of the aorta, peripheral arteries, and vena cava, respectively. The corresponding prescribed compliance of the compartments are $C_{ao}$, $C_{art}$, and $C_{vc}$.

2.1.2 Heart model

Contraction of the left ventricle is described using a time-varying elastance model [36, 37] given as

$$P_{lv}(V_{lv},t)=e(t)P_{es}(V_{lv})+(1-e(t))P_{ed}(V_{lv})\,,$$ (6)

where $P_{es}$ and $P_{ed}$ represent the end-systolic and end-diastolic pressure-volume relationship, respectively. These relationships are defined by

$$P_{es}(V_{lv}) = E_{es,\,lv}(V_{lv}-V_{lv,\,r})\,,$$ (7a)

$$P_{ed}(V_{lv}) = A_{lv}\left(\exp\left[B_{lv}(V_{lv}-V_{lv,\,r})\right]-1\right)\,,$$ (7b)

where $E_{es,\,lv}$ represents the end-systolic elastance of the left ventricle, $V_{lv,\,r}$ is the unloaded volume of the ventricle at end-systole, whereas $A_{lv}$ and $B_{lv}$ are coefficients that describe the exponential shape of the end-diastolic pressure-volume relationship (EDPVR). The function $e(t)$ describes the time course of the chamber stiffness between end systole and end diastole and is defined as follows

$$e(t)=\begin{cases}0.5(1-\cos\dfrac{\pi t}{T_{max,\,lv}});&t\leq t_{tr,\,lv}\\ 0.5\exp(-(t-t_{tr,\,lv})\dfrac{1}{\tau_{lv}});&t_{tr,\,lv}<t\leq T_{c}\,,\end{cases}$$ (8)

where $t$ is the time elapsed since the beginning of systole, $T_{max,\,lv}$ is the time-to-peak-tension, $t_{tr,lv}$ is the transition time, $\tau_{lv}$ is the relaxation time constant, and $T_{c}$ denotes the duration of a cardiac cycle.

The same relations hold for the left atrium with a time-shift, $t_{la}$, that is given by

$$t_{la}=\begin{cases}t+0.1T_{c};&t\leq 0.9T_{c}\\ t-0.9T_{c};&t>0.9T_{c}\,.\end{cases}$$ (9)

In our model, the end-systolic elastance $E_{es,\,la}$ and transition time $t_{tr,,la}$ of the left atrium are held constant, while the left ventricle’s contractility $E_{es,\,lv}$ and transient time $t_{tr,\,lv}$ are treated as model inputs. To ease the notation, we use $E_{es}$ and $t_{tr}$ (instead of $E_{es,\,lv}$, and $t_{tr,\,lv}$) to denote the contractility and transient time of the left ventricle throughout the rest of this paper.

2.2 Physics-informed neural network (PINN)

In addition to time $t$, the PINN model incorporates 10 model parameters as inputs, i.e., $R_{av}$, $R_{ao}$, $C_{ao}$, $R_{art}$, $C_{art}$, $R_{vc}$, $C_{vc}$, $R_{mv}$, $E_{es}$, and $t_{tr}$. The baseline value of these parameters are derived from previous studies [38, 39, 11]. These values are modified using a set of multipliers, which serves as the actual model inputs. Here, the allowable range for the multipliers is $[0.3,\,3.0]$ except for the multiplier of $t_{tr}$, which has a range $[0.8,\,1.2]$. Details of the baseline values and model constants are provided in Table 2 in the Appendix A. We note that while selecting a larger input range for the multipliers can enhance the model’s ability to capture a wider range of outputs, it will also complicate the training of the PINN.

Four separate neural networks, each with a size of $256\times 4$, are used here. A schematic of the network architecture is shown in Fig. 2. The output of each network represents the volume of one compartment of the closed loop system (Fig. 1). The volume of the fifth component (volume of the vena cava) is calculated based on mass conservation, i.e.,

$$\hat{V}_{vc}=V_{total}-(\hat{V}_{lv}+\hat{V}_{ao}+\hat{V}_{art}+\hat{V}_{la}).$$ (10)

We note here that PINN-based model predictions are marked with the ‘hat’ symbol (e.g., $\hat{V}_{vc}$) to distinguish them from the true values (e.g., $V_{vc}$).

To enforce the periodic boundary conditions (steady state solution), the volume waveforms are each represented by the first twelve terms of the Fourier series suggested in [40, 41]. The activation function $\tanh$ is used in all layers except the output layer.

Figure 2: Model architecture: Input layer size is 22, comprising 12 Fourier terms and 10 input parameters. Four separate neural networks, each characterized by its set of weights and biases represented as $\boldsymbol{\theta}_{i}$.

A dataset comprising $100\,000$ cases was generated for the training phase using Latin Hypercube Sampling (LHS) [42] to sample the parameter space. Each case represents a unique combination of input parameters. An additional $100\,000$ cases were selected using the same LHS method to create a test set. Each case has $400$ timepoints. In total, we have $40\,000\,000$ collocation points for training and another $40\,000\,000$ collocation points for testing.

The optimization problem for training the PINN model is defined as follows

$$\boldsymbol{\theta}^{*}=\operatorname*{arg\,min}_{\boldsymbol{\theta}}\mathcal{L}(\boldsymbol{\theta})\hskip 1.0pt,$$ (11)

where $\mathcal{L}$ is the loss function, $\boldsymbol{\theta}$ denotes $\left[\boldsymbol{\theta}_{1}^{T},\boldsymbol{\theta}_{2}^{T},\boldsymbol{\theta}_{3}^{T},\boldsymbol{\theta}_{4}^{T}\right]^{T}$, and $\boldsymbol{\theta}_{1}$, $\boldsymbol{\theta}_{2}$, $\boldsymbol{\theta}_{3}$ and $\boldsymbol{\theta}_{4}$ are the vectors of all the weights and biases for their respective networks.

The PINN model is implemented using PyTorch. The derivatives of the volume variables ($\hat{V}_{lv},\hat{V}_{ao},\hat{V}_{art},\hat{V}_{vc},\hat{V}_{la}$) with respect to time are calculated using the Automatic Derivation [43]. To solve the optimization problem (Eq. 11), the builtin ADAM optimizer was utilized with an initial learning rate of $10^{-3}$ [44]. The learning rate was dynamically adjusted during training. Specifically, if the loss stops decreasing, the learning rate was reduced by a factor of two. The minimum learning rate was set to $10^{-6}$.

To stabilize the training and prevent the model’s hyperparameters from being trapped in sharp local minima, the following strategies were used:

1. The inputs were scaled between 0 to 1 (a common practice) [45] whereas the scaling range for the outputs was adjusted to -1 and 1 that aligns with the range of the activation function $\tanh$.

2. The Mean Absolute Error (MAE) was used as the loss function to start the training. Once the learning rate becomes smaller than $10^{-4}$, the loss function was switched to the Mean Square Error (MSE). The MSE loss function is defined as the average of squared residuals for the five ODEs (Eq. 2.1.1) evaluated at all the collocation points, i.e.,

   $$\begin{split}\mathcal{L}_{MSE}(\boldsymbol{\theta})=\dfrac{1}{\mathcal{N}_{c}}\sum_{i=1}^{\mathcal{N}_{c}}&\left(\left[\frac{d\hat{V}_{lv}}{dt}-(q_{mv}-q_{av})\right]_{i}^{2}+\left[\frac{d\hat{V}_{ao}}{dt}-(q_{av}-q_{ao})\right]_{i}^{2}\right.\\ &+\left[\frac{d\hat{V}_{art}}{dt}-(q_{ao}-q_{art})\right]_{i}^{2}+\left[\frac{d\hat{V}_{vc}}{dt}-(q_{art}-q_{vc})\right]_{i}^{2}\\ &\left.+\left[\frac{d\hat{V}_{la}}{dt}-(q_{vc}-q_{mv})\right]_{i}^{2}\right),\end{split}$$ (12)

   where $\mathcal{N}_{c}$ denotes the number of collocation points. Similarly, the loss function based on MAE is defined as

   $$\begin{split}\mathcal{L}_{MAE}(\boldsymbol{\theta})=\dfrac{1}{\mathcal{N}_{c}}\sum_{i=1}^{\mathcal{N}_{c}}&\left(\left|\frac{d\hat{V}_{lv}}{dt}-(q_{mv}-q_{av})\right|_{i}+\left|\frac{d\hat{V}_{ao}}{dt}-(q_{av}-q_{ao})\right|_{i}\right.\\ &+\left|\frac{d\hat{V}_{art}}{dt}-(q_{ao}-q_{art})\right|_{i}+\left|\frac{d\hat{V}_{vc}}{dt}-(q_{art}-q_{vc})\right|_{i}\\ &\left.+\left|\frac{d\hat{V}_{la}}{dt}-(q_{vc}-q_{mv})\right|_{i}\right).\end{split}$$ (13)

   This approach produced better results than continuing the training with MAE or starting the training with MSE as the loss function because both approaches often resulted in the training algorithm being trapped in a local minima.

3. A small $\alpha$ (with a value of $0.002$) was used in the smooth maximum function (Eq. 3) to avoid the training algorithm from being trapped in the local minima at the beginning. When the learning rate becomes smaller than $10^{-4}$, the value of $\alpha$ was adjusted to get a closer approximation of the maximum function (i.e., $\alpha=0.01$).

To evaluate the performance and accuracy of the trained PINN model, a separate set of $1000$ cases was generated to compare the model results with the numerical solution of ODEs that were solved using the forward Euler Method. The evaluation metric employed here is the Relative Mean Absolute Error (RMAE), as defined by

$$\text{RMAE}(\hat{y})=\dfrac{1}{N}\sum_{i=0}^{N}\left(\dfrac{|\hat{y}-y_{N}|}{y_{N}}\right)_{i}\hskip 2.0pt,$$ (14)

where $y_{N}$ denotes the numerical solution. The model was trained on four Nvidia Tesla V100 GPUs in the High-performance Cloud Computing nodes at Michigan State University.

2.3 Sensitivity analysis

Sobol sensitivity analysis was performed to assess the influence of each model input parameter on the model outputs. The analysis decomposes the total variance in the model output into contributions from individual model input parameters and their mutual interactions [46, 47, 48, 49, 50]. To outline the mathematical framework used to calculate the Sobol indices, we denote the PINN model as a function $f$ of the input parameters $\mathbf{x}=(x_{1},x_{2},\dots,x_{s})$ that are rescaled to the unit interval $[0,1]$. Each parameter can be conceptualized as a random variable following a uniform distribution, with the assumption of mutually independence. In this probabilistic interpretation, the function $f$ is characterized as a random variable with mean $f_{0}$ and variance $V$ defined as

$$f_{0} =\int_{\Omega}f(\mathbf{x})\prod_{i=1}^{s}dx_{i}\,,$$ (15)

$$V =\int_{\Omega}f^{2}(\mathbf{x})\prod_{i=1}^{s}dx_{i}-f_{0}^{2}\,.$$ (16)

The Sobol method decomposes the total variance $V$ into summands representing contributions from individual parameters, including their interactions with each other. This is done by decomposing $f(\mathbf{x})$ into

$$f(\mathbf{x})=f_{0}+\sum_{i=1}^{s}f_{i}(x_{i})+\dots+\sum_{i=1}^{s}\sum_{j\neq i}^{s}f_{ij}(x_{i},\,x_{j})+\dots+f_{12\dots s}(\mathbf{x})\,,$$ (17)

where the decomposition terms can be calculated by following relations

$$f_{i}(x_{i})= \int_{\Omega_{\sim i}}f(\mathbf{x})\prod_{k\neq i}dx_{k}-f_{0}\,,$$ (18)

$$f_{ij}(x_{i},\,x_{j})= \int_{\Omega_{\sim(i,j)}}f(\mathbf{x})\prod_{k\neq i,j}dx_{k}-f_{i}(x_{i})-f_{j}(x_{j})-f_{0}\,,$$ (19)

By $\Omega_{\sim i}$, we denote the parameter space of all components except $x_{i}$. The higher order terms are calculated in a similar way. It can be shown that all the summands are orthogonal, which yields

$$V=\sum_{i=1}^{s}V_{i}+\sum_{i<j}V_{ij}+\sum_{i<j<k}V_{ijk}+...+V_{12\dots s}\,,$$ (20)

where

$$V_{i}=\int_{\Omega}f_{i}^{2}\prod_{k=1}^{s}dx_{k}\,.$$ (21)

The Sobol sensitivity indices are defined as

$$S_{i}=\dfrac{V_{i}}{V}\,,$$ (22)

and the total order Sobol indices are defined as

$$S_{Ti}=S_{i}+S_{ij}+...+S_{1\dots i\dots s}\,.$$ (23)

Here, only the first and second order Sobol indices were calculated, and the total order Sobol indices are the summation of the first and second order indices. The Sobol integrals were evaluated using the Monte Carlo Method with $360\,000$ functional evaluations. Time-average value of the indices, which were computed over $200$ timepoints is reported here. The SALib package in Python [51, 52] was used to compute the Sobol indices.

2.4 Inverse modeling

Once the PINN model was trained, it was then used for inverse modeling to estimate parameters from experimental measurements.

2.4.1 Experimental data

To validate the model, experimental data were acquired from 11 Yorkshire domestic swine. Left ventricular (LV) pressure waveform was acquired in each animal with a pressure catheter that was inserted through the carotid or femoral artery [53]. Left ventricular volume waveform was obtained from 3D ECHO images. These images were acquired using an EPIQ 7C ultrasound system (Philips, Andover, MA) with the animal in a supine position [54]. In 3 animals, the waveforms were acquired both before and after the administration of dobutamine, an inotropic agent that increases the myocardial contractility [55, 56, 57]. The recorded LV pressure and volume waveforms were used in the inverse modeling process to estimate model parameters and to test whether this process can quantify the increase in contractility associated with the administration of dobutamine. All animal experiments were performed in accordance with national and local ethical guidelines, including the Guide for the Care and Use of Laboratory Animals, the Public Health Service Policy on Humane Care and Use of Laboratory Animals, and the Animal Welfare Act, and an approved California Medical Innovations Institute IACUC protocol regarding the use of animals in research.

2.4.2 Parameter estimation

The goal here is to use given single-beat LV pressure and volume waveforms to estimate an optimal combination of input parameters that maximizes the coefficient of determination ($R^{2}$) between the PINN model predicted waveforms and the given waveforms. The coefficient of determination is defined by:

$$R^{2}=1-\dfrac{\sum_{i=0}^{N}(y_{true}-\hat{y})_{i}^{2}}{\sum_{i=0}^{N}(y_{true}-\bar{y})_{i}^{2}}\,,$$ (24)

where $\bar{y}$ is the average of known data, i.e.,

$$\bar{y}=\frac{1}{N}\sum_{i=0}^{N}y_{true_{i}}\,.$$

A value of 0 for $R^{2}$ indicates poor accuracy of the regression (akin to a simple mean-based model), while a value of 1 indicates that the model’s predictions capture the data perfectly. In essence, $R^{2}$ serves as a metric to evaluate how well a regression model aligns with the actual data, with higher values indicating a closer fit.

The differential evolution method [58, 59] with a population size of $100$ and a maximum of $200$ generations was used to perform the optimization. The parameter intervals employed during optimization are within the limits used for training the PINN model. The optimization was executed using SciPy [60].

To investigate the uniqueness of the model parameters, parameter estimation was also performed on $100$ pairs of synthetic waveforms generated by the PINN model with distinct sets of parameter values to determine whether the same values used for generating the waveforms are recovered with inverse modeling.

3 Results

3.1 Training of the PINN

(a) Training and test loss curves.

(b) RMAE between the PINN model and numerical results in the $1000$ evalutaion cases

Figure 3: Results of the PINN training

2(a) shows the training and test losses of the PINN model, which was trained for a total of $3300$ epochs. The noticeable decrease in losses around epoch $600$ corresponds the transition from using the MAE to using the MSE in the loss function. 2(b) shows the RMAE between the PINN model and numerical results of the 5 volume state variables (subsubsection 2.1.1) in the $1000$ evaluation cases. Average of the errors are all below $2\%$, with the largest error not exceeding $5\%$. The smallest error is found in the volume associated with the vena cava.

3.2 Sensitivity analysis

Fig. 4 shows the total order Sobol indices associated with $V_{lv}$ and $P_{lv}$ for each input parameter. The LV volume $V_{lv}$ is most sensitive to contractility $E_{es}$ followed (in order) by the vena cava compliance $C_{vc}$, peripheral artery compliance $C_{art}$, and its resistance $R_{art}$. The LV pressure waveform is also most sensitive to $E_{es}$, and is followed by $C_{art}$ and the transition time $t_{tr}$. Both $V_{lv}$ and $P_{lv}$ are insensitive to $R_{av}$, $R_{ao}$, $C_{ao}$, $R_{vc}$, and $R_{mv}$ within the parameter range applied in this study, with Sobol indices below $0.05$.

Figure 4: Total order Sobol indices associated with $V_{lv}$ and $P_{lv}$ for each input parameters

3.3 Inverse modeling on synthetic data

Fig. 5 shows the RMAE between parameters estimated based on the 100 pairs of synthetic waveforms generated by the PINN model and their ground truth. The RMAE is below $1\%$ for all parameters. We also note that the inverse model completely recovers the original waveform, where the fitted waveforms all have an $R^{2}$ exceeding $0.998$.

Figure 5: RMAE between the estimated input parameters and their ground truth for the $100$ synthetic cases

3.4 Inverse modeling on experimental measurements

Fig. 6 shows the results of inverse modeling performed using the trained PINN model with measurements from 11 swine model, including measurements acquired in 3 swine model before and after the administration of dobutamine. The fitted LV pressure and volume waveforms agree well with the experimental measurements, as reflected by the large value of $R^{2}$. Specifically, the lowest $R^{2}$ value for the LV volume waveforms and the LV pressure waveforms is $0.81$ and $0.75$, respectively.

(a) Eight swine subjects at baseline

(b) Three swine subjects before (top panel) and after (bottom) administration of dobutamine

Figure 6: Comparison of the fitted LV pressure and volume waveforms from inverse modeling with the corresponding measurements from the swine models

Table 1 shows the parameter values estimated from inverse modeling performed on the data of each swine. For the last three swine, for which measurements are available before and after dobutamine administration, the estimated $E_{es}$ (i.e., contractility) associated with dobutamine is higher. Specifically, $E_{es}$ was estimated to increase from $2.73$, $2.50$, and $1.00$ to $4.34$, $7.82$, and $3.84$, respectively, in those animals after dobutamine was administered. The transient time ($t_{tr}$) estimated by the PINN model with measurements after the administration of dobutamine is also lower compared to the values estimated with measurements before the administration.

Table 1: Estimated parameters for each set of swine data. Units for resistances, compliances, LV contractility ($E_{es}$) and transition time ($t_{tr})$ are $mmHg\>ms/ml$, $ml/mmHg$, $mmHg/ml$ and $ms$, respectively.

Animal	State	$R_{av}$	$R_{ao}$	$C_{ao}$	$R_{art}$	$C_{art}$	$R_{vc}$	$C_{vc}$	$R_{mv}$	$E_{es}$	$t_{tr}$
1	baseline	1.8	720.0	0.3	389.2	1.7	2.7	370.4	1.2	1.96	389
2	baseline	11.7	241.9	0.3	655.7	10.0	19.6	42.8	1.2	1.22	453
3	baseline	1.8	306.6	1.0	3167.0	10.0	2.7	400.0	1.2	7.09	420
4	baseline	1.8	77.4	1.0	3375.0	10.0	13.7	400.0	1.2	5.82	366
5	baseline	1.8	720.0	0.2	353.8	3.9	3.3	336.4	1.5	2.10	398
6	baseline	18.0	382.2	0.2	546.5	7.6	2.7	114.8	1.2	1.44	376
7	baseline	1.8	720.0	0.4	563.4	8.0	4.7	125.3	1.2	2.03	393
8	baseline	1.8	561.5	0.4	1862.5	1.0	5.6	400.0	4.1	3.98	567
9	baseline	5.2	436.6	0.2	817.1	1.0	7.9	181.0	12.4	2.73	345
	dobutamine	11.8	210.2	0.4	865.2	1.0	14.7	128.4	12.4	4.34	314
10	baseline	18.0	341.9	0.1	834.4	10.0	3.0	274.8	7.4	2.50	567
	dobutamine	1.8	213.9	0.1	1111.9	1.0	2.7	316.1	12.4	7.82	316
11	baseline	18.0	278.9	0.1	728.8	10.0	27.0	59.4	5.8	1.00	567
	dobutamine	1.8	322.6	0.1	1092.9	1.0	27.0	162.2	12.4	3.84	328

4 Discussion

We have developed and trained a PINN model that is based on the mathematical ODE formulation of a closed-loop lumped parameter description of the cardiovascular system. The trained PINN model can be used for inverse modeling to estimate parameters reflecting the heart function in less than 3 minutes. We also validated the inverse PINN modeling approach using experimental measurements acquired from a swine model with and without the administration of an inotropic agent.

4.1 Accuracy of the PINN model

The training strategy we adopted here enables us to minimize the loss function, resulting in an error averaging less than $2\%$ across the entire parameter space for all five volume state variables when compared to the numerical model. The largest observed error is $5\%$, which is within the level of accuracy found in other studies using PINN. Specifically, a maximum error of around $10\%$ and an average error of less than $2\%$ are reported for a PINN model developed to estimate material properties and predict the deformation of aorta [18]. Similarly, errors below $1\%$ and $7\%$ are reported for PINN models developed to simulate soft-tissue mechanics [34] and tube flows based on the Navier-Stokes equation [35]. We note here that the number of parameters used as model input in these studies (between 2 to 5 parameters) is much smaller compared to that (10 parameters) used here.

4.2 Computational time

Training of the PINN model takes 96 hours. However, the trained PINN model can be used to obtain the Sobol sensitivity analysis, which requires $360\,000$ model evaluations, within $30$ minutes. Furthermore, inverse modeling for each case, which may require up to $20\,000$ model evalutions, is completed in less than 3 minutes. We note that given that solving the ODEs numerically takes around 10 seconds to obtain the steady-state solution, performing the Sobol sensitivity analysis and the inverse modeling will require 1000 (cf. $30$ mins) and 5.56 hours (cf. $3$ mins), respectively. These results demonstrate the ability of the PINN model to efficiently perform tasks that would otherwise be significantly time-consuming if numerical methods are used.

4.3 Sensitivity analysis of the PINN model parameters

The PINN model enables us to perform a global sensitivity study efficiently where we found that $V_{lv}$ is most sensitive to $E_{es}$, followed by $C_{vc}$, $C_{art}$ and $R_{art}$ (Fig. 4). The LV pressure $P_{lv}$ is most sensitive to $E_{es}$, followed by $C_{art}$, $t_{tr}$, $C_{vc}$ and $R_{art}$. The finding that $E_{es}$, which reflects the LV contractility, affects the LV volume and pressure waveform the most is expected. The results where the compliance of the vena cava $C_{vc}$ and arteries $C_{art}$ also substantially affects the LV pressure and volume waveforms is also expected as a change in $C_{vc}$ and $C_{art}$ affects the preload and afterload imposed on the LV, respectively.

4.4 Parameter estimation using the trained PINN model

Our findings indicate that inverse modeling using the PINN model can successfully recover the generated synthetic single-beat LV volume and pressure waveforms in all the 100 evaluation cases as well as their corresponding model parameters. This result implies that these parameters are unique. Correspondingly, parameters estimated from inverse modeling using experimental measurements of the volume and pressure waveforms are likely to be also unique (Fig. 6 and Table 1), which provide confidence of using them to reflect certain physiological and pathological conditions. The smaller error in the estimation of $E_{es}$ and $t_{tr}$ parameters suggests a higher reliability in their values compared to the rest of the parameters.

The results from inverse modeling performed on the experimental measurements acquired in the animals after dobutamine administration show an increase in the estimated $E_{es}$ from $2.73$, $2.50$, and $1.00$ at baseline to $4.34$, $7.82$, and $3.84$, respectively. Correspondingly, the percentage increase of the estimated $E_{es}$ associated with the administration of dobutamine ranges between 58% to 284%. This range is higher than the values reported in previous studies, where the percentage increase is between $30\%$ to $50\%$ in Landrace pig [61], canine [62] and human [63]. This disparity may be attributed to differences in species, dobutamine dosage and the approach by which $E_{es}$ was obtained experimentally based on manipulation of the ventricular loading conditions by occluding the vena cava [61, 62] or pharmacologically using nitroprusside [63].

4.5 Potential clinical utility

The rapid inverse modeling approach described here can potentially be translated clinically to estimate ventricular contractility $E_{es}$ using measurements acquired from a single heart beat (i.e., single beat estimation). Single-beat estimation of $E_{es}$ is clinically useful because it does not require multiple pressure-volume loops with different end-systolic pressure that are usually obtained by manipulating the loading conditions that is not feasible in the clinic. While different approaches to estimate $E_{es}$ using single-beat measurements have been developed with reasonable accuracy [64, 65], they can only be used to estimate contractility but not other clinical parameters. The approach described here, which is grounded by physics, can potentially be used to simultaneously estimate other parameters (in addition to $E_{es}$) such as the systemic vascular resistance (reflected by $R_{art}$) or diastolic stiffness of the ventricle (as reflected by $A_{lv}$ and $B_{lv}$) using measurements acquired from a single beat. Further investigations (e.g., determining the type of single-beat measurements), however, are needed to realize this goal.

4.6 Limitations

This study has some limitations. First, we only consider a limited number of model parameters to reduce the computational cost associated with training the PINN model. Second, some of the estimated parameter values approached their upper or lower bounds (see Table 1), which suggests that the accuracy of inverse modeling could potentially be further improved by increasing or modifying the parameters’ range to be outside of the range at which the PINN model is trained. The parameters’ range used for training the PINN model can be expanded but this will increase the computational cost associated with the training. Nevertheless, we note that the PINN model predictions of $V_{lv}$ and $P_{lv}$ are not sensitive to those parameters that reach the bounds (e.g., $R_{mv}$ and $R_{av}$) (see Fig. 4). Last, the neural network’s architecture, which was determined in an iterative process here, may not be optimal. Determining the optimal architecture and other hyperparameters of the neural network is a notable challenge [66, 67]. Despite significant efforts in this area [68, 69], a consistent and systematic approach for selecting an optimal network architecture and determining optimal hyperparameters for PINNs remains elusive.

5 Summary and conclusions

In summary, we developed a PINN model that encodes the physics associated with the closed-loop blood circulation system consisting of the LV. The PINN model enables fast prediction of hemodynamics, which allows Sobol sensitivity analysis and inverse modeling to be executed rapidly. Parameter values of the PINN model were estimated with inverse modeling based on measurements from 11 swine models. The estimated end-systolic elastance $E_{es}$ is higher with measurements acquired after the administration of an inotropic agent in 3 swine model, which validates the inverse modeling approach. These results suggest that the PINN inverse modeling approach can potentially be further developed as a clinical tool to simultaneously estimate ventricular contractility and other physiological parameters using single-beat measurements from patients.

Code and data availability

The code and the data are available at: https://github.com/ehsanngh/lpm_pinn

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix

A Model parameters

Table 2: Constant model parameters and baseline values of the input parameters

		Parameter	Value	Unit
Baseline values for input parameters		$R_{av}$	6.0	$\mathrm{mmHg\;ms/ml}$
		$R_{ao}$	240.0
		$R_{art}$	1125.0
		$R_{vc}$	9.0
		$R_{mv}$	4.1
		$C_{ao}$	0.3	$\mathrm{ml/mmHg}$
		$C_{art}$	3.0
		$C_{vc}$	133.3
		$E_{es,lv}$	3.00	$\mathrm{mmHg/ml}$
		$t_{tr,lv}$	420	$\mathrm{ms}$
Constant Parameters	Closed-Loop	$V_{total}$	5200	$\mathrm{ml}$
		$V_{ao,r}$	100
		$V_{art,r}$	900
		$V_{vc,r}$	2800
		$T_{c}$	800	$\mathrm{ms}$
	Left-Ventricle	$A_{lv}$	1.00	$\mathrm{mmHg/ml}$
		$B_{lv}$	0.027	$1/\mathrm{ml}$
		$V_{lv,r}$	10	$\mathrm{ml}$
		$T_{max,lv}$	280	$\mathrm{ms}$
		$\tau_{lv}$	25
	Left-Atrium	$E_{es,la}$	0.45	$\mathrm{mmHg/ml}$
		$A_{la}$	0.45	$\mathrm{mmHg/ml}$
		$B_{la}$	0.05	$1/\mathrm{ml}$
		$V_{la,r}$	10	$\mathrm{ml}$
		$T_{max,la}$	150	$\mathrm{ms}$
		$\tau_{la}$	25
		$t_{tr,la}$	225

B Supplementary plots

Figure 7: First order Sobol indices associated with $V_{lv}$ and $P_{lv}$ for each input parameter

Figure 8: Second order Sobol indices associated with $V_{lv}$ and $P_{lv}$ for each input parameter

Figure 9: Total order Sobol indices associated with $V_{ao}$, $V_{art}$, $V_{vc}$, and $V_{la}$ for each input parameter

(a) PV loops for eight swine subjects at their baseline

(b) PV loops for three swine subjects before (top panel) and after (bottom) Dobutamine administration

Figure 10: Comparison of the PV loops from inverse modeling with the corresponding measurements from the swine models