Joint data imputation and mechanistic modelling for simulating heart-brain interactions in incomplete datasets

The use of mechanistic models in clinical studies is limited by the lack of multi-modal patients data representing different anatomical and physiological processes. For example, neuroimaging datasets do not provide a sufficient representation of heart features for the modeling of cardiovascular factors in brain disorders. To tackle this problem we introduce a probabilistic framework for joint cardiac data imputation and personalisation of cardiovascular mechanistic models, with application to brain studies with incomplete heart data. Our approach is based on a variational framework for the joint inference of an imputation model of cardiac information from the available features, along with a Gaussian Process emulator that can faithfully reproduce personalised cardiovascular dynamics. Experimental results on UK Biobank show that our model allows accurate imputation of missing cardiac features in datasets containing minimal heart information, e.g. systolic and diastolic blood pressures only, while jointly estimating the emulated parameters of the lumped model. This allows a novel exploration of the heart-brain joint relationship through simulation of realistic cardiac dynamics corresponding to different conditions of brain anatomy.


Introduction
Heart and brain are characterized by several common physiological and pathophysiological mechanisms [6]. The study of this multi-organ relationship is of great interest, in particular to better understand neurological diseases such as vascular dementia or Alzheimer's disease. The development of computational models simulating heart and brain dynamics is currently limited by the lack of databases containing information for both organs. Neuroimaging datasets often provide a limited number of cardiac-related measurements, usually restricted to brachial diastolic and systolic blood pressure (DBP, SBP) [7]. These quantities provide a limited assessment of the cardiac function, thus compromising the possibility of further analysis of cardiovascular factors in brain disorders.
On the contrary, the availability of a rich set of cardiac information in heart studies allows the use of cardiovascular models to estimate descriptors of the cardiac function that are not possible to measure in-vivo, such as contractility or stiffness of the heart fibers. These models optimize the parameters through dataassimilation procedures to reproduce the observed clinical measurements [11], but usually do not include neurological factors. The ability to jointly account for cardiovascular descriptors and brain information is key to gather novel insights about the relationship between heart dynamics and brain conditions.
Most of current studies relating heart and brain are based on statistical association models, such as based on multivariate regression [5]. While this kind of analysis allows to easily formulate and test association hypotheses, it usually offers a limited interpretation of the complex relationship between organs. This issue is generally addressed by mechanistic modeling, allowing deeper insights on physiological and biomechanical aspects. These models allow for example to describe the brain vasculature, and to quantify physiological aspects such as blood flow auto-regulation effects in the brain [1], up to the simulation of the whole-body circulation with detailed compartmental components [3]. Although these approaches offer a high level of interpretability, they are usually severely illposed and require large data samples and arrays of measurements to opportunely tune their parameters.
To bridge the gap between data-driven and mechanistic approaches to heartbrain analysis, in this work we propose to learn cardiovascular dynamics from brain imaging and clinical data by leveraging on large-scale datasets with missing cardiac information. This is achieved through an inference framework composed of two nested models accounting respectively for the imputation of missing cardiac information conditioned on the available cardiac and brain features, and for a Gaussian Process emulator that mimics the behavior of a lumped cardiovascular model. This setting allow us to formulate a probabilistic end-to-end generative model enabling imputation of missing measurements and estimation of cardiovascular parameters given a subset of observed heart and brain features.
Results on real data from the UK Biobank show that our framework can be used to reliably estimate and simulate cardiac function from datasets in which we have minimal cardiovascular information, such as SBP and DBP only. Moreover, the proposed framework allows novel exploration of the joint heart-brain relationship through the simulation of realistic cardiac dynamics corresponding to different scenarios of brain anatomy and damage.

Problem statement
We denote by ν the vector representing brain image-derived phenotypes (IDPs) and clinical information such as age or body surface area (BSA), and by x the vector of cardiac IDPs and blood pressure measurements. The vector x can be represented as x = {x, x obs } wherex represents the unobserved information we wish to impute and x obs the observed one. Moreover, we assume that for each observation x a corresponding set of parameters y of the associated lumped model is available. We would like to learn a generative model in which we assume that the unobserved measurements are generated by a latent random variable z conditioned on the variables ν. Hence, our generative process can be seen as sampling from a distribution p(z|ν) and then obtainingx with probability p(x|z, ν). Due to the association between cardiac IDPs, x, and lumped model parameters y, we also assume that the latter are dependent from ν through z. Our graphical model is shown in Figure 1, while the evidence lower bound (ELBO) of the joint data marginal p(y,x|x obs , ν) writes as: x y xobs Fig. 1. Graphical model of our framework. From a latent variable z we generate the unobserved featuresx conditioned on the variables ν and we estimate y via Gaussian process regression. During inference x obs and ν are used to estimate the approximated posterior q φ (z|x obs , ν).
The approximation of the posterior distribution by q φ (z|x obs , ν) defines the optimization of the ELBO through variational inference. The variational distributions q φ (z|x obs , ν) and p θ (x|ν, z) are parametrized respectively with parameters φ and θ. The term log p ω (y|x obs , z) in the ELBO denotes the emulator which approximates the mechanistic behavior of the lumped model via Gaussian Process regression parametrized by ω. The choice of a GP as emulator is motivated by the uncertainty of the data, hence it is desirable to obtain a distribution of interpolating functions rather than a single deterministic function. Moreover, GPs have already been proved to be valid emulators of 1D mechanistic vascular models [10]. The GP allow us to sample functions f (x) from a given prior parameterized by a mean µ(x) and covariance Σ(x), i.e: f (x) ∼ N (µ(x), Σ(x)) to obtain the marginal y ∼ N (µ(x), Σ(x) + σ 2 I). The prior mean µ(x) is here set to 0 and we use a radial-basis function (RBF) kernel for the covariance: Since our data is multi-dimensional, k j (x i , x ′ i ) is the kernel for the j th target and the i th predictor. The hyper-parameters of the kernel ω = {α, β} represent the output amplitude α and length scale β of the sampled functions. The goal during training is to learn the hyper-parameters that maximize the marginal likelihood of the observed data y .
The second term of the ELBO is related to the imputation ofx. The term log p θ (x|ν, z) denotes the log-likelihood of the imputed features, which can be seen as the reconstruction error. The last term KL(q φ (z|x obs , ν)||p θ (z|ν)) is the Kullback-Leibler divergence between variational approximation and prior for z, that can be expressed in a closed form given that both distributions are Gaussians. The imputation scheme is equivalent to a conditional variational autoencoder (CVAE) which has become a popular approach to feature imputation [8]. Fast and efficient optimization in our model is possible by means of stochastic gradient descent thanks to the closed form for data fit and KL terms, the use of the reparametrisation trick and Monte Carlo sampling.

Data processing and cardiovascular model
From UK Biobank we selected a subset of 3445 subjects for which T1, T2 FLAIR magnetic resonance images (MRI) and several brain and cardiac IDPs were available. Among the brain IDPs we used total grey matter (GM), total white matter (WM) and ventricles volumes. We used T1 and T2 FLAIR images to obtain the number of white matter hyper-intensities (WMHs) and their total volume relying on the lesion prediction algorithm (LPA), available from the lesion segmentation toolbox (LST) [14] of SPM 3 . WMHs are a common indicator of brain damage of presumably vascular origin [15]. We combined WM and GM volumes into a single measurement that we denoted as brain volume. The WMHs total volume and number of lesions presented a skewed distribution, and were Box-Cox transformed prior to the analysis. Regarding the cardiac IDPs we selected stroke volume (SV), ejection fraction (EF) and end-diastolic volume (EDV) for the left ventricle. All brain-related volumes were normalized by head size. Besides IDPs, we had access to blood pressure measurements (DBP, SBP) and sociodemographic features such as age and body surface area (BSA). We used DBP and SBP to compute MBP as M BP = DP + (SP − DP )/3. Next, we used the lumped cardiovascular model derived in [4] to obtain additional indicators of the cardiac function. In particular we estimated the contractility of the main systemic arteries (τ ), peripheral resistance (R p ) , the radius of the left ventricle (R 0 ), contractility of the cardiac fibers (σ 0 ) and their stiffness (C 1 ). The parameters of the model were selected based on the available clinical data. To obtain the target values for the emulator, the data-assimilation procedure was carried out according to the approach presented in [2,12].

Experiments
The data was split in two sets: one containing the full-information (2309 subjects), and one in which cardiac IDPs (x) and the estimated model parameters (y) were removed (1136 subjects, Table 1). The quality of imputation was compared to conventional methods such as mean, median and k-nearest neighbors (KNN) imputation. We subsequently assessed the relationship learnt by our model between brain features and cardiovascular parameters through simulation. Starting from the mean values of each parameter we sampled along the dimension of the different conditional variables ν that we used to parameterize the prior. This procedure allowed us to assess their influence in the inferred simulation parameters. Table 1. Variables used in our framework. Mean blood pressure (MBP), diastolic blood pressure (DBP), stroke volume (SV), end-diastolic volume (EDV), ejection fraction (EF), heart fibers contractility (σ0), ventricle size (R0), heart fibers stiffness (c1), peripheral resistance (Rp) and aortic compliance (τ )

Results
Data imputation and regression. We assessed the performance of our model by measuring the mean squared error (MSE) on the testing data, for both the imputation of missing cardiac information and the emulation of the lumped model parameters. The results in Figure 2a show that our method gives significantly better estimates than mean and median imputation, and comparable predictions to KNN, for which the optimal number of neighbors was optimized through 10fold cross-validation and corresponds to K = 10. At the same time the emulator consistently gives low errors for the parameters' estimation ( Figure 2b). In supp. Figure 3 we present the most relevant predictors for each emulated feature based on their β values, while a qualitative comparison between the distributions of imputed and emulated features compared to the ground truth data is available in in supp. Cardiovascular dynamics simulation and model plausibility. In Figure  3a we observe the change in the predicted parameters of the cardiovascular model as we sample along the range of values of the different conditional variables ν. Figure 3b shows the pressure-volume (PV) loops generated by the lumped model using the inferred parameters. The simulated PV loops highlight meaningful relationships: -An increase in the volume of WMHs is associated to decreased SV and EDV, together with a smaller reduction of ESV, leading to a decrease of EF which is related to reduced contractility. -Similar dynamics are associated to brain volume loss.
-The number of WMHs exhibits different dynamics than WMHs volume, associated to the increase in afterload and the increase of EDV.
From the plots showing the estimation of the parameters we can observe the ones driving the observed dynamics. For example, WMHs and brain volumes changes are mainly driven by the joint evolution of peripheral resistance (R p ), contractility (σ 0 ) and size of the left ventricle (R 0 ). The changes in the number of WMHs are related to heart-remodelling changes, driven by R 0 and by a decrease of σ 0 and the ventricular stiffness c 1 . In supp. Figure 2 we provide analysis and discussion of the remaining conditional features (age, BSA and brain ventricles volume). Overall, in the simulated dynamics we can identify several physiological responses in line with the clinical literature. Moreover, our model may be useful to give insights in currently controversial topics such as concerning the pathogenesis of WMHs. Our results suggest that the evolution of cardiac function with respect to brain and WMHs volumes is similar to the one due to aging (see supp. similar to the one related to ventricles enlargement. In both cases the changes in R 0 , σ 0 and C 1 suggest that the increase in the number of WMHs and the enlargement of the brain ventricles are related to heart-remodelling processes, such as loss of contractility or decrease in compliance. These findings are in line with clinical observations [9] relating lower cardiac output with higher burden of WMHs and reduced brain volume.

Conclusion
We presented a generative model that enables the analysis of complex physiological relationships between heart and brain in datasets where we have minimal available features. The framework allow us to emulate a lumped cardiovascular model through data-driven inference of mechanistic parameters, and provides us a generative model to explore hypothetical scenarios of heart and brain relationships. In the future, the model will allow to potentially transfer the knowledge learnt in UK Biobank to datasets where we have minimal cardiac information, to explore the relationship between brain conditions and cardiovascular factors in specific clinical contexts, such as in neurodegeneration. Our approach could also be extended to account for deep learning architectures, and the framework could be further improved by jointly accounting for multiple outputs, which are currently modelled independently, or by including spatial information from imaging data, beyond the modelling of scalar volumetric features. Furthermore, while the cardiovascular features considered in this study are rather general, more complex features will allow to study more realistic cardiovascular models. For example, while the mechanistic model used in this study does not simulate cerebral blood flow, previous studies suggested that WMHs may be due to local vascular impairment [13]. Hence, by selecting appropriate clinical features we could constrain the imputation by means of any biophysical model representing the desired aspect of systems biology. This could represent an innovative tool in real world scenarios, for which multi-modal patient data is often limited or not available.