Classification of myocardial blood flow based on dynamic contrast‐enhanced magnetic resonance imaging using hierarchical Bayesian models

Dynamic contrast‐enhanced magnetic resonance imaging (DCE‐MRI) is a promising approach to assess microvascular blood flow (perfusion) within the myocardium, and the Fermi microvascular perfusion model is widely applied to extract estimates of the myocardial blood flow (MBF) from DCE‐MRI data sets. The classification of myocardial tissues into normal (healthy) and hypoperfused (lesion) regions provides new opportunities for the diagnosis of coronary heart disease and for advancing our understanding of the aetiology of this highly prevalent disease. In the present paper, the Fermi model is combined with a hierarchical Bayesian model (HBM) and a Markov random fields prior to automate this classification. The proposed model exploits spatial context information to smooth the MBF estimates while sharpening the edges between lesions and healthy tissues. The model parameters are approximately sampled from the posterior distribution with Markov chain Monte Carlo (MCMC), and we demonstrate that this enables robust classification of myocardial tissue elements based on estimated MBF, along with sound uncertainty quantification. A well‐established traditional method, based on a Gaussian mixture model (GMM) trained with the expectation–maximisation algorithm, is used as a benchmark for comparison.

estimated MBF, along with sound uncertainty quantification. A well-established traditional method, based on a Gaussian mixture model (GMM) trained with the expectation-maximisation algorithm, is used as a benchmark for comparison.

K E Y W O R D S
classification, dynamic contrast-enhanced magnetic resonance imaging, hierarchical Bayesian model, Markov random fields prior, myocardial blood flow, myocardial perfusion INTRODUCTION Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is a non-invasive method for the assessment of myocardial perfusion (see Jerosch-Herold, 2010). This method is based on the intravenous administration of MRI contrast agent, which selectively changes the signal intensity in proportion to the local contrast agent concentration in the myocardium. In this way, DCE-MRI can capture signal intensity changes in the time domain during the first pass of the contrast agent through the myocardium.
Several methods for the analysis of myocardial perfusion DCE-MRI scans have been proposed in the literature (e.g. Jerosch-Herold, 2010). Most of these methods aim to quantify myocardial blood flow (MBF), but do not address the problem of classification of myocardial tissue with respect to perfusion level and, hence, disease status, which is indispensable for the accurate diagnosis of coronary heart disease (CHD). The ultimate objective of the present work is to pave the way towards a clinical decision support system for automatic lesion detection, which is based on quantifying the difference between the MBF of normal (healthy) and hypoperfused (lesion) tissues. To this end, we use the Fermi microvascular perfusion model (see Jerosch-Herold et al., 1998), which has been widely applied in the literature to estimate the MBF. However, we emphasise that this particular choice of MBF model is not critical, and that it is straightforward to combine the modelling framework we propose with other tracer-kinetic MBF quantification methods, including tracer-kinetic models (see Biglands et al., 2015;Larsson et al., 1996).
Our work is based on hierarchical Bayesian modelling, which has been extensively reviewed, for instance, in chapters 5 and 15 of Gelman et al. (2014). Hierarchical Bayesian models (HBMs) have been applied in a broad and diverse range of research and development areas, including healthcare medicine (see Kantorová et al., 2020), transportation networks (see Pitombeira-Neto et al., 2020), clinical trials (see Chu & Yuan, 2018), water resources (see Bracken et al., 2018) and economics (see Meager, 2019). Moreover, HBMs have been applied widely in medical imaging, for example DCE-MRI. A HBM method has been developed to estimate kinetic parameters of DCE-MRI based on Gaussian Markov random fields (MRF) priors in Schmid et al. (2006Schmid et al. ( , 2009. introduced a HBM method to model the concentration-time curve (CTC) based on a standard compartment model (see Parker, 2005) where a traditional nonlinear regression model was used for comparison. Whitcher et al. (2011) illustrated an application of Schmid et al. (2009) to evaluate treatment response of squamous cell carcinoma. Of particular interest in the context of the present article are recent applications of HBMs to myocardial perfusion modelling. Schmid et al. (2007) modelled the non-linear attenuation of arterial input functions (AIFs) and used HBMs to reconstruct AIFs from magnetic resonance (MR) signal intensities in the myocardium. Schmid (2011) extended the modelling in Schmid et al. (2007) by importing spatio-temporal constrains to the B-spline model using Gaussian MRF. Lehnert et al. (2019) illustrated an HBM method to analyse the myocardial perfusion DCE-MRI. In this method, a large-scale Bayesian spatio-temporal regression method was applied to model the data. Similarly, Gaussian MRF priors were introduced to express the spatio-temporal constrains. Scannell et al. (2020) used HBMs to infer the parameters of tracer-kinetic models from contrast-enhanced MR images. They further combined HBMs with MRF to include prior knowledge on the physiological ranges of the model parameters and to integrate additional information on the physiological neighbourhood of a target voxel. The commonality of all these applications is that the models adopt hierarchical structures to incorporate potentially complex interactions and dependencies between their variables, and use state of the art Markov chain Monte Carlo (MCMC) sampling techniques to address the analytical intractability of inference.
Our work aims to provide an automated identification of areas of hypoperfused myocardium. To this end, we combine a myocardial perfusion model with a novel HBM and a MRF prior. The latter acts on discrete labels that we associate with the pixels of an MRI scan, to assign them to two distinct clusters: lesion versus healthy tissue. Loosely speaking, this clustering process can be viewed as a two-layered model, where the first (observed) layer represents the signal intensity (proportional to contrast agent concentration) and the second (hidden) layer contains the class labels (see Sui et al., 2012). A MRF prior is used to provide context information, so that class assignment is not only based on a pixel's own signal intensity, but also takes spatial information from neighbouring pixels into account. As we will demonstrate, this is essential for noise reduction, good generalisation and robust classification, that is obtaining classification patterns that are unaffected by random noise injections.
We note that the application of HBMs to myocardial perfusion modelling has been proposed in four previous publications: Schmid et al. (2007), Schmid (2011), Lehnert et al. (2019) and Scannell et al. (2020). Here, the emphasis has been on the estimation of the MBF per se, and the authors have demonstrated that the MBF parameters can be more accurately estimated than with traditional methods, such as non-hierarchical non-linear regression, particularly in regimes where the signal-to-noise ratio (SNR) is low (see Scannell et al., 2020 for details). The essential novel contribution of our work is the additional introduction of discrete labels associated with the individual image pixels, and their regularisation via a MRF. This methodological extension is essential for automatic tissue classification and lesion detection.

BACKGROUND
In DCE-MRI, the tissue of interest (e.g. myocardium) is imaged before, during and after intravenous administration of exogenous MRI contrast agent (see Figure 1). In T1-weighted DCE-MRI, gadolinium-based contrast agents (GBCA) selectively shorten local T1 (spin-lattice or longitudinal) relaxation times, in proportion to local tissue concentration of GBCA (C(t)). This gives rise to changes in post-injection DCE-MRI signal intensities (S(t)) compared to baseline (S(0)): relatively hypoperfused regions remain dark during the first-pass of GBCA, compared to regions characterised by higher perfusion (see Figure 1). Local tissue concentration of the contrast agent at time t, C(t), can be approximated by the difference between post-contrast signal S(t) and baseline signal S(0): where t denotes time. Any departure from this simplified assumption of linearity will affect the absolute values of extracted parameters, but will not affect their relative distribution as long as S(t) is a monotonously rising function of C(t) over the relevant range of concentrations within the myocardial tissue. With the contrast dose used in this study (0.1 mmol∕kg), MBF values will be systematically overestimated. This is why we refer to them as 'estimates of MBF' rather than 'absolute MBF' within the context of this report.

Quantification of myocardial blood flow with the Fermi method
Methods for estimation of MBF can be separated into two categories: model-based and model-independent. In this section, we will briefly review the latter approach, of which the Fermi model used in our work is a particular example. For a more comprehensive review, including the former category, the reader is referred to Jerosch-Herold (2010). The model-independent methods are based on the central volume principle, which was introduced in Zierler (1965). This principle can be mathematically described as where F b is the MBF and C myo (t), C in (t), C out (t) denote the contrast agent concentrations in the myocardium, arterial input (concentrations in arterial blood) and venous output (concentrations in venous blood) at time t, respectively. For any linear and stationary system, C out (t) can be expressed by a convolution of C in (t) and a transfer function h(t) (see Zierler, 1965 for details), in which the symbol * denotes the convolution operation. The transfer function h(t) denotes the distribution of the transit time required for the contrast agent molecules to traverse the myocardial tissue from the arterial blood pool. Substituting Equations (3) into (2) leads to: where R f (t) denotes the impulse response function, whose expression is given by The details of the derivation of Equation (4) can be found in Jerosch-Herold et al. (1998).
Once R f (t) is deconvolved based on Equation (4), F b can be obtained by setting t = 0 in R f (t) as indicated in Equation (5). In general, a unique and universal form of R f (t) does not exist, and its free-form inference is an 'ill-posed' problem due to its intrinsic unidentifiability (see Jerosch-Herold, 2010). In order to estimate F b , usually a function family is chosen for R f (t), such as the widely used Fermi function (see Axel, 1983;Jerosch-Herold et al., 1998), piece-wise smooth B-splines (see Jerosch-Herold et al., 2002) subject to Tikhonov regularisation (see Tikhonov et al., 2013), exponential deconvolution (see Keeling et al., 2009) and autoregressive moving average (ARMA) models (see Neyran et al., 2002;Batchelor et al., 2010). Since these models have parameters without a direct physical interpretation, they are typically referred to as 'model-free' (which in the context of statistical modelling is a misnomer). This is to be distinguished from model-based approaches, whose parameters do have a direct physical interpretation. In the present work, we use the Fermi model, which is a 'model-free' approach whose representation of the impulse response function has the following parametric form: where A, , are shape parameters. We obtain the MBF from the amplitude of R f (t) at t = 0. Our clinical data sets are limited to first-pass of the contrast agent through the vasculature, with very restricted coverage of the extravasation. This is why we chose to implement the Fermi model, as opposed to compartmental modelling approaches (see Biglands et al., 2015). We emphasise again, though, that the method we propose can work, in principle, with any model of R f (t), and the choice of the Fermi model is not critical to our approach. Traditionally, least-squares fitting has been used to estimate the parameters of MBF models (see e.g. Jerosch-Herold et al., 1998;Scannell et al., 2020). This method aims to minimise the cost function, where t denotes temporal index of the dynamic image, i is the pixel index and is the vector of pixel-dependent Fermi parameters. C(t| i ) is the model-predicted contrast concentration and y i (t) is the observed contrast concentration, both at pixel i and time t. However, the algorithm of this optimisation method may get stuck in local minima (see Dikaios et al., 2017;Scannell et al., 2020). Moreover, merely finding the parameters that minimise Equation (7) fails to quantify the uncertainty intrinsic to the parameter estimation problem. For that reason, we follow a Bayesian approach and aim to approximately sample the parameters from the posterior distribution using MCMC sampling (see e.g. Gelman et al., 2014). In this way, we avoid entrapment in local minima and properly quantify the uncertainty inherent in the parameter estimation. We note that a similar approach has been adopted in Scannell et al. (2020). The details of our MCMC sampling scheme are provided in the next section.

The application of Markov random fields prior
As mentioned before, our objective is the DCE-MRI pixel classification for automatic lesion detection. However, the naive application of standard clustering methods to the DCE-MRI pixel intensities can lead to spurious singleton clusters, especially when the SNR in the data is low. This also applies to simple generative models, like the Gaussian mixture model, as we will demonstrate in the Results section. Such singleton clusters, or more generally small fragmented clusters, are physiologically unrealistic due to the pathology of myocardial ischemia (see Yang et al., 2019). To address this problem, we use a MRF prior. This prior uses spatial context information based on the assumption that adjacent pixels reflect similar cardiophysiological kinetic properties. We will show in the Results section that this will lead to physiologically more realistic MBF estimation and myocardium pixel classification.
Let k i ∈ {0, 1} denote the label of pixel i. When k i = 0, this pixel is marked as 'healthy tissue'. When k i = 1, this pixel is marked as 'lesion'. According to the Hammersley-Clifford theorem, a MRF is equivalent to a Gibbs distribution (see Li, 2009), where Q is a normalisation constant to ensure that P(k i |k −i ) is normalised and sums to 1, k −i is the set of all other labels except k i , k −i = {k j } j≠i , and in which T k i is a weight parameter, i ∼ m means pixel i and pixel m are neighbouring pixels and u(k i |k m ) is a local potential function, which we define to be where 'o' indicates the neighbouring degree of order. In this paper, we will use up to second-order neighbourhoods, o ∈ {1, 2}. The definition of spatial pixel cliques (pixel i and its neighbours) can be found in Figure 3. Increasing the neighbourhood of a pixel in the MRF to third order was found to lead to a substantial increase in the computational complexity, by about 100%. The restriction to second order is motivated by the fact that in nearly all our images, the thickness of the myocardium wall does not exceed five pixels. A third-order neighbourhood could potentially straddle the myocardium wall and combine regions on either side of it, which is an undesirable effect. The main justification for our restriction to second order, however, is our empirical finding that increasing the neighbourhood from second to third order was not found to make any difference in practice. Figure 2 shows a comparison between segmentations obtained with a second-order order and a third-order neighbourhood, leading to very similar results. where

For the Fermi parameters
and Note that the dependence of the right-hand side of Equation (13) on k i and k −i is via Equation (15). In Equation (13), T Fermi is a weight parameter, and i ∼ m indicates that pixels i and m are neighbouring pixels, which are defined in Figure 3. The parameter family { i,m } is used to guarantee edge-preserving properties that is to prevent smoothing over edges that are indicative of tissue boundaries (e.g. delineating the area of a lesion); see Scannell et al. (2020) and Bardsley (2012) for details. When two adjacent pixels are in different categories, the parameters { i,m } give the model the flexibility to switch off the interaction between the pixels and ensure they are not going to influence each other in the MBF. To be specific, in this work, i,m is defined by

Hierarchical Bayesian model
As we have discussed in the previous section, we use binary labels to assign pixels to two classes, related to 'healthy tissue' and 'lesion'. According to the pathology of myocardial perfusion, healthy tissue can be expected to have substantially higher MBF than lesions. We can model this with a prior distribution of the MBF parameters in Equation (6) where the parameters of this distribution, the so-called hyperparameters, are dependent on the binary labels. We apply the hierarchical Bayesian modelling framework (Gelman et al., 2014) to consistently model the hierarchical dependence among all parameters, hyperparameters and labels. The structure of this HBM is a directed acyclic graph (DAG). The DAG used in the present work is shown in Figures 4 and 5. A DAG allows the following factorisation of the joint probability distribution (see e.g. Bishop, 2006):

F I G U R E 3
where K is the total number of nodes and pa q is the set of parents of q . The parents are the nodes with an edge (arrow) feeding into the target node. For example, for parameter i in Figure 4, the parents nodes are , −i , k i and k −i . This expansion rule is used in Section 3.3 to factorise the joint distribution of the HBM. In Figures 4 and 5, if there is a link from node A to node B, then node A is the parent of node B, and the child node is conditionally dependent on the parent node. The descriptions of the parameters in Figures 4 and 5 can be found in Table 1. The structure of the HBM model reflects the three-layer causal relationships of the data-generating process. Specifically, the observed data is dependent on the contrast agent and its concentration, as described by the Fermi model, and the effect of the contrast agent is dependent on the status label. The rationality of the former dependence is from the central volume principle (see e.g. Zierler, 1965) and the Fermi model (see e.g. Jerosch-Herold et al., 1998), and the rationality of the latter dependence follows from F I G U R E 4 This figure shows the structure of the hierarchical Bayesian model proposed in the present article, where we assume that the impulse response function is a Fermi function. The parameter groups are defined below Equation (21). In this graph, the higher layers are conditionally dependent on the lower layers. The circle nodes denote variables and the rectangle nodes denote fixed values or observations. Figure 4 in which all parameters are shown individually. This inevitably leads to a more cluttered diagram, which has the advantage, though, that the conditional probabilities required for the Gibbs sampling scheme of Section 3.4 can directly be read off the graph. Recall that the conditional probability of a target node conditional on the complete set of random variables in the graph is equal to the conditional probability of the target node conditional on its Markov blanket. The Markov blanket is the set of parents, co-parents and children, where a parent is a node with a directed edge pointing to the target node, a child is the recipient of a directed edge from the target node, and the co-parents are the other parents of the target node's children. For more details, see chapter 8 in Bishop (2006).

Parameter
Description

Means of the prior distributions of the Fermi parameters for pixel
Variances of the prior distributions of the Fermi parameters for pixel i Labels of the neighbours of pixel i * and ** The symbols with * and ** are hyperparameters the physiological fact that the MBFs of healthy tissues and lesions are different (see e.g. Knaapen et al., 2006). The details of the nodes in Figures 4 and 5 will be illustrated in Section 3.3.

Novel statistical model
Let y i (t) denote the observed contrast agent concentration of an image where i is a positive integer (i = 1, 2, … , N) and t is a non-negative integer (t = 0, 1, … , M − 1). N is the total number of pixels in this image, M is the total number of time stamps of the data. For each pixel, k i denotes the group of this pixel. In this case, we assume that there are two groups, healthy tissue and lesion. Therefore, k i = 0 when the ith pixel is in the healthy group and k i = 1 when the ith pixel is in the lesion group. We assume the observed contrast concentration y i (t) is Gaussian distributed with mean C(t| i ) and variance 2 i . Therefore, the expression of y i (t) is where C(t| i ) is the model calculated contrast agent concentration. Note that the notation C(t| i ) means the same as C myo (t) in Equation (4), but we have made the dependence on the Fermi parameter vector i explicit in the conditioning argument. i ∼ N(0, 2 i ) is assumed to be i.i.d. additive Gaussian noise (see Chatzis & Varvarigou, 2008;Zhang et al., 2001). We will test in our simulation studies how critical this distributional assumption is for the ultimate purpose of this study (automatic lesion detection). The parameters i are the Fermi parameters defined in Equation (8). According to the central volume principle (see Zierler, 1965, section 2.1), the expression of contrast agent concentration C(t| i ) can be found in Equation (4), which we rewrite to make the dependence on the Fermi parameters, defined in Equation (8), explicit in our notation: where R f (t, i ) can be found in Equation (6). The new modified notation makes the dependence on the Fermi parameter vector i explicit in the argument. The term C in (t) in Equation (18) can be directly acquired from data, as explained in Section A of the online supplementary material. We now make the dependence on the model parameters explicit in our notation. Since we assume the noise i to be iid Gaussian distributed, we have We assume that the parameters i are independent of time t. The quantities i are related to the patient's physiology and pathophysiology, which will change with time (e.g. as a consequence of treatment). However, these changes are slow, typically happening in the order of weeks or months. Our MRI time series, on the other hand, are acquired within a minute, that is, several orders of magnitude below. It is therefore justified to assume that for the duration of the measurement modelled in our study, the parameters i can be taken as constant. For As we have discussed in Section 3.2, according to the factorisation rule, the joint distributions of a DAG is the product of conditional distributions determined by the respective parent nodes. Given the structure in Figure 4, the joint distribution is The joint probability distribution in Equation (21) has combined two classes of probabilistic graphical models: directed probabilistic graphical models, also called Bayesian networks (reviewed e.g. in Murphy, 2012, section 10), and undirected graphical models, also called MRF (reviewed e.g. in Murphy, 2012, section 19). Both modelling paradigms can be integrated into the generalised concept of a partially directed acyclic graphical model, also called chain event graph or just chain graph; see e.g. Koller and Friedman (2009) or the recent monograph by Collazo et al. (2018). This combination of a Bayesian network with a MRF is the natural approach for the present problem. The DAG structure of the Bayesian network represents the causal relationships of the data-generating process, as explained in the previous section. Spatial neighbourhood relations, on the other hand, lack a natural DAG structure, and an undirected graph, that is a MRF, is a more natural representation. Our model can thus be interpreted as a particular instance of a chain graph, one in which the directed and undirected graph structures are clearly separated, and this separation simplifies the inference compared to the more general inference schemes discussed in Collazo et al. (2018).
We now explain each term on the right-hand side of Equation (21). The likelihood P(y i | i , 2 i ) can be found in Equation (20).
is the prior distribution of the Fermi parameter vector i , which combines explicit priors for the three parameter groups with the MRF prior of Equation (12): ) are Gaussian prior distributions. Pathologically, the values of MBF for healthy tissues and lesions are different. Since all Fermi parameters are non-negative, we assume the logarithms of the Fermi parameters are Gaussian distributed conditional on the labels k i . We have also tried Gamma distributions as an alternative, and include the results in Section B of the online supplementary material. The difference in the results obtained with a log-Gaussian versus a Gamma prior was found to be minor, which suggests that the choice of functional family for the prior distributions on the Fermi parameters i is not critical, as long as the distributions are consistent with the positivity constraint of the Fermi parameters, that is have positive support.
we have: with k i ∈ {0, 1}. For i and i , the definitions are similar. (12) and (13). We choose a vague inverse gamma distribution for the observational noise variance 2 i , due to its conjugacy. P( | * * , k i ) is the prior distribution for hyperparameter vector . It can be written as where and with ∈ {0, 1}. { * * = 0, 2 * * = 10, * * = 0.1, * * = 0.1} are uninformative hyperparameters. P(k i |k −i ), another term in Equation (21), is the MRF prior of label k i , which can be found in Equations (9)-(11). P( −i ) and P(k −i ) in Equation (21) are the marginal distributions of −i and k −i . While these expressions are intractable, they do not enter the Gibbs sampling scheme for posterior inference, to be discussed in the next section.

Posterior inference
We now derive a Gibbs sampling scheme for posterior inference. In a DAG (the HBM model in Figure 5), the probability of a selected target parameter conditional on all other parameters is given by the probability of the target parameter conditional on its Markov blanket (see e.g. chapter 8 in Bishop, 2006). The Markov blanket is the set of parents, children and co-parents. For example, for the parameter 2 i in Figure 4, its parent nodes are * and * , its child node is y i (t) and its co-parent is i . In this way, the relevant conditional probabilities can immediately be read off the graph in Figure 5 in Section 3.2.
For the parameter vector i = {A i , i , i }, according to Figure 4, its parents nodes are −i , , k i and k −i . Its child node is y i (t) and its co-parent node is 2 i . Therefore, its total conditional posterior distribution is: where P(y i | i , 2 i ) can be found in Equation (20) and P( i | −i , , k i , k −i ) can be found in Equation (22). Note that P(y i | i , 2 i ) depends on C(t| i ); see Equation (20). C(t| i ) is defined by a convolution integral-Equation (18)-which is analytically intractable. Hence, the conditional posterior distribution of i in Equation (28) is analytically intractable.
According to Figure 4, the parent, child and co-parent nodes of 2 i are * , * , y i (t) and i . The conditional distribution of the observational noise variance 2 i , is obtained by combining the likelihood for pixel i, with the conjugate prior for P( 2 i | * , * ) from Equation (24) to give where * = * = 0.1, from the previous section. According to Figure 4, the parent, children and co-parent nodes of k i are k −i , i , , * * and −i . Its conditional posterior distribution is given by can be found in Equation (22); P( | * * , k i ) can be found in Equation (25), and P(k i |k −i ) is the MRF prior for binary parameters, and its expression can be found in Equations (9)-(11). The conditional posterior distribution of k i is a Bernoulli distribution because k i ∈ {0, 1}. Specifically, when k i = 0, we have ) .
Note that the term i,m depends on the label k i given Equation (15). Similarly, when k i = 1, we have ) .
In this way, the conditional posterior distribution for label k i is According to Figure 5, the parent nodes, child node and co-parents nodes of A,k i and 2 A,k i are { * * , 2 * * }, { * * , * * }, k i , A i , A −i and k −i . Therefore, the conditional posterior distributions are given by and where ∈ {0, 1} is the binary lesion indicator. As seen from Equation (22), the term that includes the dependence on the neighbourhood, A −1 and k −i , does not depend on A,k i and 2 A,k i . Therefore the terms A −i and k −i drop out. By simplifying Equation (36), the conditional posterior distribution for A,k i = is where A,k i = . By simplifying Equation (37), the conditional posterior distribution for 2 A,k i = is The conditional posterior distributions for { ,k i = , 2 ,k i = , ,k i = , 2 ,k i = } can be derived similarly.
In the Bayesian statistics framework, samples from the conditional posterior distributions of i , 2 i , k i , and can be approximately obtained with MCMC methods. In this work, the conditional posterior distribution of i is analytically intractable because of the convolution form in Equation (18). Therefore, we adopt a Metropolis-Hastings-within-Gibbs sampling scheme to sample the parameters. To be specific, i is sampled with the Metropolis-Hastings algorithm based on adapted Gaussian proposals distributions, as discussed in Chapters 11 and 12 of Gelman et al. (2014). All other parameters will be sampled directly from their conditional posterior distributions by Gibbs sampling. The details of our sampling method can be found in algorithm 1.

Clinical data
Four patients (2 male, age range 56-60), underwent a comprehensive cardiac MRI exam on a Siemens MAGNETOM Avanto (Erlangen, Germany) 1.5-Tesla scanner with a 12-element phased array cardiac surface coil. DCE-MRI assessment of resting myocardial perfusion was performed during intravenous administration of 0.075 mmol∕kg of contrast agent (gadoterate meglumine, Dotarem, Guebert S.A.). Myocardial tissue segmentation (definition of endocardial and epicardial contours as well as extraction of LV blood pool data) was performed using cardiac image analysis software (QMASS 8.1; Leiden, The Netherlands).
In the following analysis, we choose one patient to illustrate our method (56-year-old male, scan acquired 2.7 days after percutaneous coronary intervention (PCI), with left anterior descending (LAD) artery myocardial infarction). The results obtained in the remaining three data sets are presented in the online supplementary material.

Synthetic data
Since it is challenging to obtain the ground truth of the size and location of the lesion in clinical data, we also use synthetic data to test our method. Two assumptions are made. First, the noise for data generation is assumed to be either Gaussian or Rician noise. Due to the central limit theorem, Gaussian noise is the most general and widely applied form of noise, while Rician noise is often assumed for MRI data due to the underlying physical characteristics (see e.g. Gudbjartsson & Patz, 1995). Note that by generating data subject to Rician noise, we test the robustness of our method (which assumes Gaussian noise, as discussed in Section 3.3) with respect to model mis-specification. Second, the change of signal intensity for each pixel is a double exponential function, as explained below; see Equation (40) for the temporal variation of indicator concentration in the peripheral compartment (see Radjenovic, 2003), which mimics the change of signal intensity of each pixel. The synthetic data are designed to emulate the clinical data we have used in this work: 46 images (the same number as for the clinical data) are simulated to mimic the acquired DCE-MRI data. The first 12 images are pre-contrast images, and the remaining 34 images are contrast enhanced images. We discard the first 12 images. The first image of the remaining 34 images corresponds to time point t = 0. We assume that there is one lesion, which is located in the anteroseptal region of the myocardium. The other pixels in the myocardium mimic healthy tissue with normal blood perfusion. For the pre-contrast image, the signal intensities for all pixels are drawn from a Gaussian distribution : N(10, 1). These values have been chosen to ensure that, at this stage, all signal intensities are designed to be low. For the contrast enhanced images, the modelling of the signal intensities is more complex. We use a double exponential curve to mimic the change of signal intensity with time, which is given by the following equation: We use different parameters for healthy tissue and lesions. For lesion pixels, the parameters are set to p 1 = 0.02, p 2 = 0.3, p 3 = 20, while for healthy tissue pixels, the parameters are set to p 1 = 0.01, p 2 = 0.4, p 3 = 25. These values are chosen empirically to mimic our clinical data. The process to obtain these values can be found in Figure 6. For Gaussian noise, the signal intensity in the contrast enhancement stage is s(t) + N(10, 2 * ). The variance of the noise 2 * is chosen to be 2 or 2.5, corresponding to SNRs of 4.1 and 3.72, respectively. For Rician noise, we use the Rayleigh distribution to generate the noise; see Gudbjartsson and Patz (1995) for details. The value 2 * is chosen as 1.8 with an SNR of 2.91. Again, all these values have been selected so as to mimic our clinical data.

Alternative models for comparison
We have compared the performance of the HBM model proposed in the present paper with two alternative models for estimating the MBF. The first alternative model predicts signal intensities independently according to the Fermi model, as defined in Equation (6), without taking any neighbourhood information of the image pixels into account. The parameters are estimated by minimising the residual sum-of-squares score between measured and predicted signals, defined in Equation (7). We refer to this model as the Fermi model. The second alternative model is a simplified HBM model where no explicit discrete labels for lesion indication are included. This model is similar to the model proposed in Scannell et al. (2020). The only difference is that we have replaced the authors' MBF model with the Fermi model of Equation (6), for the reasons already discussed in the text below Equation (6). We refer to this model as the HBM model without labels.
The model parameters are sampled from the posterior distribution with MCMC, using the same sampling scheme as described in Scannell et al. (2020). We have also included a comparison with two alternative classification methods, for assigning pixels to the classes 'healthy tissue' versus 'lesion'. The first alternative scheme is a standard Gaussian mixture model, whose parameters are estimated in a maximum likelihood F I G U R E 6 The process of generating realistic synthetic data based on real clinical data. See Section 4.2 for more detailed explanations. [Colour figure can be viewed at wileyonlinelibrary.com] sense using the expectation maximisation (EM) algorithm. Building a standard Bayes classifier, this leads to a well-established classification method that is widely used in the literature (see e.g. Bishop, 2006). We apply this method to the MBF map predicted by the Fermi model, described in the previous paragraph, and refer to the combined scheme as the GMM based on Fermi method. The second alternative classification method is an extension of the method just described with a morphological noise removal correction called the 'opening and closing operation' (see Chen & Haralick, 1995 for details). We refer to this method as GMM + O&C.

Model selection and convergence diagnostics
The parameters of our model are sampled from the posterior distribution using the sampling algorithm described in Section 3.4. The exceptions are hyperparameters T k i and T Fermi defined in Equations (10) and (13). These hyperparameters define the strength of the neighbourhood influence in the MRF. They could in principle be sampled from the posterior distribution using MCMC methods, along with the other parameters. However, this would cause a considerable increase in the computational complexity, because convergence and mixing of lower-level hyperparameters tends to be happening at much lower rates than for higher-level parameters (see e.g. Neal, 1996 for more details). For clinical decision making, excessive computational costs need to be avoided. We have therefore chosen discrete values for these parameters, with equidistant intervals on a log 10 scale. These values were obtained from pilot studies where the influence of the neighbourhood information in the MRF was explored. We then treat the selection of these parameters as a model selection problem, using the Watanabe information criterion (WAIC) introduced in Watanabe (2013); see also chapter 7 in Gelman et al. (2014) for more textbook-style explanations. This score, which can be directly computed from the MCMC trajectories, is asymptotically equivalent to Bayesian leave-one-out cross-validation, with lower values indicating models that are preferred by the data. The advantage is that we can now run MCMC simulations for different values of T k i and T Fermi in parallel, and then select the parameters with the lowest WAIC scores. This leads to substantial computational savings compared to directly sampling these parameters from the posterior distribution with MCMC.
To assess the convergence of the MCMC simulations described in Section 3.4, we computed the Geweke statistic (see e.g. Geweke, 1991 ) directly from the MCMC trajectories, taking Z-scores in the range of (−1, 1) as absence of evidence for insufficient convergence.

Results for synthetic data
We randomly generated 10 independent instantiations of the synthetic data, described in Section 4.2, to increase the reliability of our performance assessment. Figure 7 shows an example of the estimated MBF for the synthetic data. For each data set, we calculated the relative errors of the MBF estimations, E MBF , as follows: where N total is the total number of pixels, F i est is the estimated MBF for the ith pixel and F i truth is the ground truth of the MBF for the ith pixel. Lower values of E MBF are indicative of better performance. A comparison of the proposed HBM model with the two alternative methods for MBF estimation described in Section 4.3 is shown in Table 2. This table confirms quantitatively the qualitative impression already obtained from Figure 7: our method clearly outperforms the Fermi model, and performs slightly better than the HBM without labels. The reason for the performance improvement over the Fermi model is the inclusion of spatial context information via the MRF prior. The slight performance improvement over the HBM without labels demonstrates the edge-preserving properties afforded by the labels.    Figure 8 shows an example of automated tissue classification (lesion vs. healthy tissue). We computed the marginal posterior probability of the pixel labels from the MCMC trajectories aŝ where M 0 is the length of the Markov chain and k p i is the pth sample of label k i . Recall that k i = 1 indicates that the ith pixel is in a lesion, whereas k i = 0 indicates that the pixel lies in a healthy region. We also apply a decision threshold of 0.5 for automatic binary pixel classification. We compare the classification obtained with the new proposed HBM method with, first, the ground truth and, second, with the two alternative pixel classification schemes summarised in the final paragraph of Section 4.3. The figure suggests that the classification obtained with the proposed HBM is in very good agreement with the ground truth. The proposed HBM also clearly outperforms the GMM based on Fermi method. This is mainly due to the MRF prior, which achieves a substantial noise reduction. The opening and closing method, GMM + O&C, discussed in Section 4.3, also achieves a considerable noise reduction (see panel c in Figure 8). However, this noise reduction is not as reliable as the one achieved with our HBM model (panel d). In particular, the opening and closing method still leads to a spurious lesion, which is avoided by the proposed HBM. To follow the visual inspection up with quantitative analysis, we did two things, both exploiting the fact that we have the ground truth labels. First, we used the estimated posterior probabilities of the pixel class labels to assign each pixel to one of the two classes (0 vs. 1), based on a decision threshold of 0.5 (i.e. pixels with posterior probability above 0.5 are assigned to class 1, otherwise they're assigned to class 0). This crisp assignment does not take posterior uncertainty into consideration, i.e. the posterior probabilities 0.51 and 0.99 are both mapped to the same class. To account for the different degrees of uncertainty, we therefore also computed the cross-entropy between the true labels and the predicted posterior probabilities. The results are shown in Tables 3 (misclassification) and 4 (cross-entropy), and clearly demonstrate that the proposed HBM outperforms the two alternative methods on both scores.

Results for clinical data
In this section, we discuss the application of the proposed method to the clinical data described in Section 4.1. We start with a discussion of our MCMC convergence diagnostics, and then continue with model selection, MBF visualisation and automated pixel classification.

MCMC convergence
Recall that our MCMC sampling algorithm is described in Algorithm 1. Figure 9 shows MBF traceplots for four typical pixels. The MCMC acceptance ratios were controlled within the range (0.23, 0.44), which has been suggested to give the best mixing and convergence (see Gelman et al., 2014). The traceplots show the exploration of the MBF samples for each pixel. The traceplots for the MBF estimations of the four typical pixels look reasonable, with no trends to indicate insufficient convergence. As a more objective and quantitative convergence measure, we have computed the Geweke's statistic (see Geweke, 1991). The results are shown in Figure 10, based on four independent MCMC simulations. Since the Geweke's statistic stays consistently within the uncritical range (−1, 1), there is no significant evidence for lack of convergence. Figure 9 also shows the marginal posterior distribution of the MBF, which we obtained from the MCMC trajectories with a standard kernel density estimator. This provides a natural way to quantify the MBF estimation uncertainty inherent in the data.

Model selection
The MRF prior, introduced in Section 3.1, depends on two hyperparameters, T k i and T Fermi ; see Equations (10) and (13). Rather than sampling them from the posterior distribution, we adopted a model selection approach: we ran the MCMC simulations for various fixed values of T k i and T Fermi , and selected the best combination with WAIC; see Section 4.4 for details. The results are shown in Table 5, from which the best parameter combination can be obtained.  The bold number indicates the smallest WAIC value which represents the best combination of the hyperparameters.

Marginal posterior distribution analysis
The posterior distributions are particularly interesting in comparison with the corresponding prior distributions. For illustration, we have randomly selected a pixel (i = 399) in the MR image, and estimated the posterior distributions for the associated parameters A i , A i and 2 i . The reason we choose these three parameters is that they are representative of the other parameters. i and i are related to A i , as these are all Fermi parameters. 2 A i is related to A i ; both are hyperparameters of the prior distribution of A i . Panels (a)-(c), Figure 11 show a comparison of the corresponding prior and posterior distributions. From the strong collapse of the posterior, whose effective support has shrunk considerably compared to that of the prior, we learn that the data are highly informative for estimating these parameters. Model developers interested in interpreting the model parameters can learn about their degree of practical identifiability from these plots.
Panel (d), Figure 11 shows violin plots of the marginal posterior distributions of Fermi parameters A i , i , i and variance 2 i given a randomly selected pixel i = 399. These plots allow a quantification of the intrinsic posterior estimation uncertainty. In addition, they can be inspected for various patterns in the posterior distributions. For instance, the marginal posterior distributions of the parameters A i , i and i are symmetric and bell-shaped, consistent with a Gaussian distribution. The marginal posterior distribution of 2 i , on the other hand, is skewed with a long upper tail. This is consistent with the prior distribution (log normal or gamma distribution) and suggests that despite the fact that the data are highly informative about this parameter, with a strong collapse of the posterior distribution, as discussed above, the functional form of the prior is maintained. 4.6.5 Automated pixel classification Figure 12 shows estimations of the MBF. This still requires manual inspection and interpretation. What we need for clinical decision support is automated classification of pixels and regions as healthy versus diseased, so as to automatically detect lesions. To this end, we compute the posterior probabilities of the labels from the MCMC trajectories using Equation (42). Hence, a high posterior probability close to 1 is indicative of healthy tissue, whereas a value close to 0 indicates a lesion. Figure 13 shows the posterior probability of the class label for different models and different settings of the hyperparameters T k i and T Fermi . The figure shows a crisp separation into regions with a high posterior probability of being a lesion (shown in dark green), and regions with a high posterior probability of being healthy (shown in yellow). There is only a small 'uncertain' region with posterior probabilities in the order of 0.5, which we have enlarged in panel (f). For automatic pixel classification, we have to threshold these posterior probabilities: classify a pixel i as lesion ifP(k i = 0|D) ≥ 1 − ; classify a pixel i as healthy tissue ifP(k i = 0|D) ≤ ; and classify a pixel i as uncertain if <P(k i = 0|D) < 1 − . To investigate the effect of this threshold, we have chosen the hyperparameters with the best WAIC score (T k i = 0.1 and T Fermi = 1) and repeated the  Figure 14, which demonstrates that the classification labels are not affected by the variation of the threshold. This is in line with Figure 13, which shows that most pixels have posterior probabilities close to 1 or 0 and are therefore not affected by the particular choice of threshold. Figure 13 shows the posterior probabilities for both the GMM based on the traditional Fermi model (as described in Section 4.3), and the proposed HBM model with various settings of the hyperparameters T k i and T Fermi . Note that the traditional method in Figure 13 panel (a) shows a fragmented structure, with various small and singleton lesions embedded in healthy tissue. This is physiologically unrealistic. Left ventricular myocardium is supplied by three major epicardial arteries, and any impairment in blood delivery is expected to affect consecutive regions of tissue (vascular territories), and not isolated small patches as shown in Figure 13(a).
Interestingly, the physiologically unrealistic fragmentation and the prediction of spurious singleton clusters are suppressed with our HBM model, irrespective of the hyperparameter settings. This is shown in Figure 13 panels (b-e). This finding does not only demonstrate that the inclusion of spatial context information via the MRF prior achieves physiologically more realistic cluster assignments, but also demonstrates that the classification results show reasonable robustness with respect to a variation of the hyperparameters. Different values of the hyperparameters do not affect the location of the detected lesion much. Again, the best classification map can be found with our model selection scheme based on WAIC, shown in Table 5, which corresponds to panel Figure 13  We have run an additional simulation study to focus on the interactions between the parameter and label uncertainties. We have chosen three pixels in distinct characteristic positions of the myocardium: the centre of the healthy tissue, the centre of the lesion, and the edge between healthy tissue and lesion. We set the Fermi parameters to their 10% and 90% posterior quantiles, as obtained from our previous MCMC simulations. All other parameters (including the Fermi parameters associated with the other pixels) were kept fixed. We ran 1 MCMC simulation with 1000 Gibbs sampling steps given different quantiles of Fermi parameters (specific pixel). We selected two combinations of Fermi parameters that represent the lowest MBF and the highest MBF respectively. The results can be found in Figure 15. It can be seen that the labels of the pixel inside the healthy tissue (marked by a blue circle) and the lesion (marked by a green circle) do not change and are invariant with respect to changing the quantiles of the Fermi parameters. However, the pixel near the boundary between the healthy tissue and the lesion (marked by a yellow circle), is affected by the quantiles of the Fermi parameters, with different quantiles leading to different labels. To explain this behaviour, note that the effect of the Markov random field prior on an individual target pixel is strong when the surrounding pixels have the same label, leading to low transition probabilities into different labels for the target pixel. A target pixel near a boundary, on the other hand, is surrounded by pixels with different labels, whose influence effectively cancels out in the Markov random field. Consequently, the Fermi parameters now have a significant effect on the target pixel's label. given fixed parameters (10% and 90% quantiles of the Fermi parameters) in different locations (healthy region, lesion region, boundary region) of the myocardium successively. The white pixels indicate healthy tissue and the red pixels indicate lesions. The pixels inside the blue circles (healthy tissue), green circles (lesion) and yellow circles (boundary) have fixed Fermi parameters (10% and 90% quantiles). [Colour figure can be viewed at wileyonlinelibrary.com]

DISCUSSION AND FUTURE WORK
The principal bottleneck of the proposed method for clinical applications is the high computational cost of the MCMC simulations. In our work, a single MCMC iteration took in the order of 2 s, leading to run times of about 5.5 h for a typical MCMC chain length of 10,000 iterations (using the following hardware: Intel(R) Core(TM) i9-7900X CPU @ 3.30 GHz processor with 64 GB memory). In our future work, we will explore the potential of parallel processing to reduce these run times. This will exploit the fact that all relevant computations are invariant with respect to a permutation of the order in which the pixel is processed, and all pixel-related operations can be carried out in parallel rather than sequentially. A consequent improvement over CPU-based simulations would be the porting of our code to GPU platforms. In this way, the expected run time could be reduced to 1∕10 according to Sui et al. (2012). While MCMC is intrinsically sequential, there are alternative sampling methods that can better exploit the scope for parallelisation, most notably sequential Monte Carlo (SMC) and adaptive multiple importance sampling (AMIS), which in combination with GPU architectures could significantly reduce the run times (see Cornuet et al., 2012 andDoucet et al., 2001). Another approach is to replace sampling based Bayesian inference with modal or distributional approximations, using the Laplace approximation, variational inference or expectation propagation; see chapter 13 in Gelman et al. (2014) for details. A further limitation is related to the physical model used for quantifying the MBF (discussed in Section 2.1). The MBF is approximated by the Fermi model, which depends on three shape parameters {A, , }. For modelling the dependence on other kinetic parameters, like fractional plasma volume and fractional interstitial volume, more complex models need to be employed. However, as we have already mentioned before, the hierarchical Bayesian modelling framework we have proposed is agnostic to the specific form of the physical model employed. To paraphrase this: if a more complex physical model is employed, it can just take the place of the Fermi model in the proposed hierarchical Bayesian model, and all its parameters can be estimated methodologically the same way as proposed here for the Fermi model (exploiting the information-sharing capacity and uncertainty quantification of the hierarchical Bayesian model and its spatial Markov random field prior). To be more specific, if a model satisfies the central volume principle in Equation (2) and its structure of impulse response function R f (t) is given, this model can be put into the hierarchical structure presented in this work by replacing the Fermi parameters A, , in Figure 5 with the parameters in its impulse response function R f (t). In particular, the classification of the myocardial tissue for automatic detection of lesions will proceed in the same way as described here, while potentially benefitting from the richer structure of more advanced physical models.
Finally, there is scope for future improvement related to the choice of the prior distribution. In the present work, we have used vague prior distributions for the hyperparameters for Fermi parameters, that is the parameters in Figure 4. This is a consequence of an intrinsic limitation of the Fermi model: while widely used in the literature, it depends on empirical shape parameters that lack any clear physical interpretation. However, for more advanced and realistic physical models, for example the tracer-kinetic models (see Sourbron & Buckley, 2011a, 2011b, whose parameters have bio-physical and physiological interpretability, methods of cardio-physiological knowledge elicitation can be adopted to derive more informative prior distributions that have the potential to boost the predictive and diagnostic performance of the proposed methodology.

CONCLUSION
In this work, we have introduced a hierarchical Bayesian model with a MRF prior to classify myocardial tissues based on the MBF. This method systematically combines signal intensities from raw MR images with spatial context information related to the individual image pixels to achieve physiologically more realistic myocardial tissue classification. Traditionally, least-squares fitting has been (and is still widely being) used to estimate MBF using the Fermi model, and we have used it as a benchmark for the MBF estimation in the present work. Our work is methodologically related to, and has been inspired by, the HBM model of myocardial perfusion MRI proposed in Scannell's work Scannell et al. (2020). However, this study does not address the problem of automatic myocardial tissue classification, and we have found in general that this topic has only been addressed scarcely in the current cardio-physiological modelling literature. For that reason, we have used a standard and widely applied statistical model as a further benchmark for tissue classification from MBF data: the Gaussian mixture model, with parameter estimation using the EM algorithm. A comparison with these two traditional methods has revealed three advantages of the proposed hierarchical Bayesian modelling framework. First, it can generate a clearer, smoother and more realistic MBF map by taking physiological context information into account. Second, our inference method properly quantifies the uncertainty in both the MBF as well as the latent label assignment on which the myocardial tissue classification is based. Thirdly, the method exploits the relationship between the MBF estimation and the latent pixel labels for automatic myocardial tissue classification. This achieves automation (which is essential for clinical applications) at improved estimation and classification accuracy, paving the way to a future clinical decision support system.