Physics-informed neural networks for myocardial perfusion MRI quantification.

Tracer-kinetic models allow for the quantification of kinetic parameters such as blood flow from dynamic contrast-enhanced magnetic resonance (MR) images. Fitting the observed data with multi-compartment exchange models is desirable, as they are physiologically plausible and resolve directly for blood flow and microvascular function. However, the reliability of model fitting is limited by the low signal-to-noise ratio, temporal resolution, and acquisition length. This may result in inaccurate parameter estimates. This study introduces physics-informed neural networks (PINNs) as a means to perform myocardial perfusion MR quantification, which provides a versatile scheme for the inference of kinetic parameters. These neural networks can be trained to fit the observed perfusion MR data while respecting the underlying physical conservation laws described by a multi-compartment exchange model. Here, we provide a framework for the implementation of PINNs in myocardial perfusion MR. The approach is validated both in silico and in vivo. In the in silico study, an overall decrease in mean-squared error with the ground-truth parameters was observed compared to a standard non-linear least squares fitting approach. The in vivo study demonstrates that the method produces parameter values comparable to those previously found in literature, as well as providing parameter maps which match the clinical diagnosis of patients.

Introduction

The evaluation of dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) has proven to be a popular tool for the assessment of tissue physiology in various diseases (Nagel, Klein, Paetsch, Hettwer, Schnackenburg, Wegscheider, Fleck, 2003, Turnbull, 2009, Heye, Culling, Valdés Hernández, Thrippleton, Wardlaw, 2014). In particular, stress perfusion cardiac magnetic resonance (CMR) is becoming an established technique for the non-invasive assessment of patients with suspected coronary artery disease (CAD) (Montalescot et al., 2013). A series of large trials have demonstrated its efficacy in clinical practice (Greenwood, Maredia, Younger, Brown, Nixon, Everett, Bijsterveld, Ridgway, Radjenovic, Dickinson, Ball, Plein, 2012, Schwitter, Wacker, Wilke, Al-Saadi, Sauer, Huettle, Schönberg, Luchner, Strohm, Ahlstrom, Dill, Hoebel, Simor, 2013, Greenwood, Ripley, Berry, McCann, Plein, Bucciarelli-Ducci, Armellina, Prasad, Bijsterveld, Foley, Mangion, Sculpher, Walker, Everett, Cairns, Sharples, Brown, 2016), and it has recently been shown to be non-inferior to the invasive reference standard, fractional flow reserve (FFR), for the management of patients with suspected CAD (Nagel et al., 2019) and to be cost-effective (Kwong et al., 2019). However, since the images are difficult to interpret and ischaemic burden may be underestimated for patients suffering from multi-vessel disease, visual assessment of the severity of CAD using stress perfusion CMR remains limited (Villa et al., 2018). As such a quantitative analysis of perfusion is proposed as a reproducible and user-independent alternative to the visual assessment (Patel et al., 2010).

The quantification of myocardial blood flow (MBF) has shown promising diagnostic accuracy and prognostic value (Lockie, Ishida, Perera, Chiribiri, De Silva, Kozerke, Marber, Nagel, Rezavi, Redwood, Plein, 2011, Hsu, Jacobs, Benovoy, Ta, Conn, Winkler, Greve, Chen, Shanbhag, Bandettini, Arai, 2018, Sammut, Villa, Di Giovine, Dancy, Bosio, Gibbs, Jeyabraba, Schwenke, Williams, Marber, Alfakih, Ismail, Razavi, Chiribiri, 2018) and has the potential to allow more widespread clinical adoption of stress perfusion CMR. The basis for myocardial perfusion quantification is the use of tracer-kinetic modelling, which allows the relation of the concentration-time curves derived from DCE-MRI and patient-specific physiological parameters, such as MBF (Sourbron and Buckley, 2012). In the case of myocardial perfusion this is achieved by modelling how the contrast agent passes from the left ventricle (LV) into the myocardium, allowing for the inference of kinetic parameters.

The identification of the kinetic parameters is not, however, trivial. Several studies have shown that reliable quantification of myocardial perfusion was only achieved with relatively simple models such as Fermi-constrained deconvolution, but not with multi-compartment exchange models (Broadbent, Biglands, Larghat, Sourbron, Radjenovic, Greenwood, Plein, Buckley, 2013, Schwab, Ingrisch, Marcus, Bamberg, Hildebrandt, Adrion, Gliemi, Nikolaou, Reiser, Theisen, 2015, Likhite, Suksaranjit, Adluru, Wilson, Di Bella, 2017). One of the underlying reasons for the difficulties in model fitting is that non-linear regression problems tend to get stuck in local optima (Kelm, Menze, Nix, Zechmann, Hamprecht, 2009, Dikaios, Atkinson, Tudisca, Purpura, Forster, Ahmed, Beale, Emberton, Punwani, 2017). As a result, the model-based concentration curves may match the noisy observed data, but the inferred parameters can be far from the actual values. It has further been shown that model parameters are correlated (Romain et al., 2017), which causes several distinct parameter combinations to produce concentration curves which may all well fit the observed data, therefore making the identification of the correct set of parameters difficult (Buckley, 2002, Ahearn, Staff, Redpath, Semple, 2005). This has lead to more complex fitting algorithms being proposed but these have not yet seen widespread adoption (Kelm, Menze, Nix, Zechmann, Hamprecht, 2009, Dikaios, Atkinson, Tudisca, Purpura, Forster, Ahmed, Beale, Emberton, Punwani, 2017, Scannell, Chiribiri, Villa, Breeuwer, Lee, 2020, Dikaios, 2020).

Given the limitations of the currently used methods, this study introduces a novel class of algorithms for the quantification of myocardial perfusion: physics-informed neural networks (PINNs) (Raissi et al., 2019). PINNs are based on the universal approximation theorem of neural networks which is leveraged to solve supervised learning tasks while respecting the given physical laws in terms of ordinary or partial differential equations (Hornik et al., 1989). Specifically, the solutions to the tracer-kinetic model differential equations are approximated by neural networks which are trained to produce outputs that both fit the available data, and satisfy the underlying physical conservation laws.

In this work, we investigate the use of physics-informed neural networks for the estimation of kinetic parameters from stress perfusion CMR data in both simulated data and a small patient cohort.

Methods

Tracer-kinetic modelling

Tracer-kinetic models present a mathematical description for the physics of the transport process of a contrast agent across a tissue (Ingrisch and Sourbron, 2013). The two-compartment exchange model (2CXM) has been suggested as appropriate for modelling myocardial perfusion (Jerosch-Herold, 2010). In this case, the perfusion unit (represented by a single pixel in DCE-MRI) is modelled as a system of two interacting compartments, the plasma and the interstitium. This gives rise to a pair of coupled ordinary differential equations (ODEs) which describe the evolution of the concentration of contrast agent over time:Here, and are the concentration of contrast agent, gadolinium [Gd] in the plasma and interstitial space at time t, respectively (in units of molarity (m)). , the arterial input function (AIF), is the assumed input to the system that is being modelled (also in m). In myocardial perfusion quantification this is sampled from the LV. is the plasma flow (ml/min/ml), is the fractional plasma volume (dimensionless), is the fractional interstitial volume (dimensionless) and PS is the permeability-surface area product (ml/min/ml). A weighted sum of the concentration within these separate compartments then produces a representation for the concentration within the myocardial tissue ( in m):This is fit to the observed imaging data to infer the parameters: , and .

Physics-informed neural networks

PINNs are a new framework within the deep learning paradigm which employ deep neural networks to approximate the solution to physical systems, primarily ordinary and partial differential equations (Raissi et al., 2019). With this framework, the parameters of the neural network are constrained and learned in two ways. Firstly, the direct outputs of the neural network are trained to fit the observed data, and secondly, the neural network is constrained to satisfy the underlying physical laws that govern this observed data. In this particular case, the physical laws are modelled by the ODEs derived from the 2CXM, found in Eqs. (1) and (2).

As such, a neural network is defined that approximates the solution to these equations with parameters for all pixels in an imaging slice. This is the mapping from time to the solutions of the ODEs at time :The outputs of the neural network are also constrained by an additional loss function derived from the set of residuals corresponding to the 2CXM:As well as learning the parameters of the neural network, these residuals may then also be used to learn the kinetic parameters present in the ODEs. This ODE residual loss encourages the PINN to produce physically plausible results and is made possible by the automatic differentiation functionality of deep learning frameworks. A schematic representation of the PINN is shown in Fig. 1.

Fig. 1

A schematic representation of the physics-informed neural network showing the process used to infer both the neural network model parameters and the tracer-kinetic model parameters . For a given time point , the neural network approximates , and . The approximations of and are combined with to approximate as in Eq. (3). These approximations (the green nodes) are compared to observed MRI concentrations and this difference is the loss term . The residuals of the ODEs are then computed (the orange nodes) and summed to give . Then, and (along with any additional regularisation terms) are minimised with respect to (w.r.t) and to find the optimal parameters.

The loss function for training the neural network is then given as the sum over all pixels of the weighted sum of these ODE residuals (), the sum of squared differences between the predicted and the observed concentrations (), the initial conditions () which enforce that there is not contrast in the system at time ), and a regularisation term to enforce non-negative concentration values ():, , , and are the weightings for each term in the combined loss, is the total number of pixels, and the exact form of each term in Eq. (6) is given in Appendix A.

Reduced form 2CXM

In this work, the idea of constraining the residuals of the total myocardial concentration rather than the plasma and interstitial compartment concentrations (as in Eqs. (1) and (2)) will be tested. This can be achieved by taking the time derivative of Eq.  (3) and replacing its terms with the ODEs of the 2CXM found in Eqs. (1) and (2):This ODE is further rewritten into a residual function, that can be incorporated into the PINN:

Experiments

The proposed algorithm will be tested in two separate ways on both simulated and patient data, allowing for a qualitative and quantitative assessment regarding the benefits of the methodology compared to the currently used non-linear least squares (NLLS) solution.

Generation of the digital phantom

In the first case a 2CXM digital reference object (DRO) is constructed, for which the pixelwise parameters are known. 144 myocardial tissue concentration curves with different parameter combinations of , , , and are subsequently simulated using a gamma-variate function for the AIF, producing the data which the model should fit. This is based on the DRO proposed by Debus et al. (2019) with the tissue kinetics, AIF, and time resolution adapted to more realistically represent myocardial perfusion. The 2CXM parameters used for this study are therefore as follows:This DRO gives rise to a volume of time curves with spatial dimensions and a unique subset of 2CXM parameters for every block of pixels. This is simulated at a temporal resolution of 0.02 min over a time span of 2 min, creating tissue curves with a total of = 100 time points as a result. Random normal noise was subsequently added to the generated curves such that the final curve has a signal-to-noise ratio of 17.5.

Digital phantom study

Using the DRO, a comparative study of four separate methods is performed, each with the goal of estimating the 2CXM parameters given the generated noisy data. The different versions of the 2CXM PINN correspond to the optimisation of different residual terms in the loss function.

The 2CXM PINN in which the original two residual functions provided in Eq. (5) are minimised in conjunction with the described PINN framework.

The Reduced 2CXM PINN which solely relies on minimising the reduced-form ODE residual found in Eq. (8).

The Combined PINN which jointly optimises the 2CXM and the Reduced 2CXM residuals. Since the residuals given both in Eqs. (5) and (8) should both be minimised, it follows that the sum of these residuals should also be minimised. Therefore the residual loss term of the Combined PINN is the sum of the residuals from the 2CXM PINN (Eq. (5)) and the Reduced 2CXM PINN (Eq. (8)).

The standard NLLS optimisation method.

For each of these methods, the normalised mean square error (NMSE) between the true and estimated kinetic parameters is reported along with a structural similarity index (SSIM) to test the overall structural coherence as compared to the ground truth (GT). The metrics were computed for each of the three 2D images and the study was repeated for 5 different noise realisations with the mean (standard deviation) performance metrics reported. The described performance metrics are used to compare the different PINN-based methods with the best performing method being used to test the feasibility of PINN-based kinetic parameter estimation in vivo, in the next section.

Patient data study

The proposed method was tested on clinical stress perfusion scans for 8 patients. The study was conducted in accordance with the Declaration of Helsinki (2000) and was approved by the National Research Ethics Service (15/NS/0030). All patients provided written informed consent. The patient datasets used for this project comprise of examinations performed on a 3T system (Achieva TX, Philips Healthcare, Best, The Netherlands), with a 32-channel cardiac phased array receiver coil. For each patient, a total of three LV short-axis slices were attained at the apical, mid-cavity, and basal level. These slices were acquired at mid-expiration with a saturation-recovery gradient echo method. Stress images were obtained during adenosine-induced hyperaemia, and for each acquisition 0.075 mmol/kg of bodyweight gadolinium contrast agent was administered at 4 mL/s, followed by 20 mL saline flush per acquisition. Each bolus of gadobutrol was preceded by a diluted pre-bolus with 10% of the dose in order to mitigate the non-linear relationship between MR intensity values and contrast agent concentration (Ishida et al., 2011).

All acquired perfusion images were subsequently corrected for respiratory motion using a previously described image registration scheme (Scannell et al., 2019b). The regions of interest for the AIF and myocardium were segmented automatically using an existing deep learning-based pipeline (Scannell et al., 2020b). The AIF is averaged over a region of interest chosen in the LV blood pool and pixelwise signal-intensity curves were extracted from the myocardial segmentation, and were split into time intervals corresponding to the pre-bolus injection and the main bolus injection in order to perform quantification. Due to the use of a dual-bolus acquisition, a linear relationship can be assumed between the signal intensity values and contrast agent concentrations. Thus, the conversion to concentration values is approximated using the relative signal enhancement approach (Ingrisch and Sourbron, 2013). Finally, the American heart association (AHA) representation was used to assign the pixelwise parameter estimates to AHA segments through automatically computed right ventricular insertion points (Cerqueira et al., 2002). The haematocrit value (HCT) was assumed to be 0.45 and the specific density of the myocardium was assumed to be 1.05 g ml. These were used to convert plasma flow and volume ( and ) to blood flow and volume ( and ).

The kinetic parameters of these patients are estimated using the best model found in the above Digital phantom study. Qualitative assessment of the stress perfusion scans was also conducted by an experienced Level-3 cardiologist (A.C) (Plein et al., 2011). The distribution of the MBF values in AHA segments which were visually positive for ischaemia was compared to the MBF values in segments negative for ischaemia by means of a boxplot. The NMSE between the AIF predicted by the neural network and the patient’s measured AIF is also reported.

Optimisation details

The traditional non-linear least squares fitting was implemented using the SciPy implementation (Virtanen et al., 2020) of the L-BFGS nonlinear optimisation scheme to minimise the residuals between the model and the observed concentration data. The optimisation was constrained to ensure all kinetic parameters are positive, and and are less than 1. For the PINN parameter fitting, in order to enforce the estimation of positive values for the physiological parameters, the logarithm of the parameters are optimised during training and the exponential of these values returned as the parameter estimates. These parameters are further initialised to a value range which is within the physiological range. The loss function is also adapted to penalise predictions of negative concentration values. While the loss on the observed measurements is enforced at the available observed measurements, the residual solution is enforced at a total of = 500 points generated through a random uniform distribution within the time domain. The ability to enforce the residuals at these collocation points is the reason why the AIF is also predicted by the neural network, as in Eq. (4).

The neural network weights are initialised using Glorot uniform initialisation and are iteratively updated using the Adam optimisation algorithm with a learning rate of 0.01, decreased by a factor of 10 every 10,000 epochs (Glorot, Bengio, 2010, Kingma, Ba, 2015). The network uses fully-connected layers and consists of two hidden layers with 32 units each, which are all followed by a hyperbolic tangent (tanh) activation function and a batch normalisation layer. This choice of activation function is justified by the choice of relatively shallow neural networks commonly employed, and the importance of balanced gradient flow required (LeCun et al., 2012). Models are trained for 25,000 iterations. The architecture was inspired by previous work in the field (Raissi et al., 2019) and there was no attempt to optimise this. The number of iterations was chosen empirically along with the weight factors , , , and . All aforementioned methods and optimisations are implemented with the Tensorflow 2 deep learning library (Abadi et al., 2016).

As is common for deep learning, the input and output data are normalised. The input data is standardised as follows:where and are the mean and standard deviation of the temporal coordinates respectively. Secondly, the output concentration values are normalised to the range in order to further mitigate any scaling performed by the neural network itself:The ODEs are re-written to account for these scalings, as shown in the Appendix B. The use of normalised inputs and outputs to the neural network and batch normalisation layers are known to protect against vanishing (or exploding) gradients and to stabilise the training process (Glorot and Bengio, 2010).

The open source implementation is provided on github1, as well as scripts to reproduce the results. All derived data is provided, including the trained model weights, but the raw patient data is not shared for ethical reasons.

Results

A summary of the DRO study performance is presented in Table 1, showing both the mean (standard deviation) NMSE and SSIM between the estimated parameters and the ground-truth DRO. These metrics are shown for each of the individual kinetic parameters, as well as the total value averaged over all the four kinetic parameters. An example 2D slice inference of the 3D DRO is presented in Fig. 2. This figure demonstrates a slice taken from the z-axis, along which the PS value varies, therefore producing a figure for which the GT value of the PS parameter is constant (). It is seen that there is a lower overall NMSE for both the 2CXM PINN and the Combined PINN as compared to the standard NLLS approach. Indeed, the Combined PINN is the best performing model overall and for all kinetic parameters except for . This model is therefore chosen for use in the patient data study results presented in Section 4.2.

Table 1

The NMSE for the individual and combined kinetic parameter estimates, evaluated for each proposed methodology. The SSIM value for the estimates is presented in brackets next to each NMSE score. Bold values represent those scores which were the best amongst all models. Scores are based on the digital reference dataset with SNR of 17.5 described in the Methods section.

	Total		F_p		v_p		v_e		PS
	NMSE	SSIM	NMSE	SSIM	NMSE	SSIM	NMSE	SSIM	NMSE	SSIM
2CXM	0.13 (0.09)	0.53 (0.15)	0.19 (0.09)	0.75 (0.19)	0.04 (0.02)	0.47 (0.11)	0.01 (0.01)	0.58 (0.12)	0.30 (0.25)	0.34 (0.16)
Reduced	0.20 (0.10)	0.24 (0.04)	0.08 (0.10)	0.73 (0.06)	0.04 (0.004)	0.20 (0.08)	0.17 (0.05)	0.03 (0.01)	0.51 (0.29)	0.01 (0.01)
Combined	0.11 (0.09)	0.57 (0.11)	0.12 (0.06)	0.80 (0.03)	0.03 (0.02)	0.49 (0.14)	0.01 (0.01)	0.62 (0.12)	0.26 (0.26)	0.37 (0.14)
NLLS	0.17 (0.18)	0.50 (0.14)	0.03 (0.02)	0.62 (0.03)	0.07 (0.07)	0.34 (0.15)	0.02 (0.02)	0.43 (0.25)	0.55 (0.61)	0.50 (0.14)

Fig. 2

Results for the inference of the tracer-kinetic parameters taken from a single 2D slice of the DRO for a single noise realisation, performed by the four methods described in Section 3.1. Inference maps are provided for all four parameters, with the first columns denoting the ground truth maps. The next four columns present the results for the different PINN methodology, with the last row showing the results for the analytical solution. NMSE and SSIM values are provided for each map relative to the ground truth.

5/8 patients were reported as being visually positive for ischaemia and 3/8 patients were reported as having visually normal images. Table 2 shows the median value (25th percentile, 75th percentile) for all inferred kinetic parameters. These results are broadly in line with values reported in the literature previously (Scannell et al., 2020a).The estimated MBF maps for two representative patients (one with and one without ischaemia) are shown in Fig. 3 along with the MR images. The corresponding images for all other patients are shown in Appendix C, along with a representative comparison of the predicted AIF with the measured AIF. Homogenous blood flow is shown in the normal case with areas of clearly reduced blood flow seen in the ischaemic patient. Fig. 4 shows the evolution of the constituent loss terms and the mean value (over the whole patient) of the estimated kinetic parameters over the training process for a representative image slice (patient 3, basal slice). The MBF () median value (25th percentile, 75th percentile in AHA segments with ischaemia on the visual assessment was 1.11 (0.82, 1.45) ml/min/g. The equivalent value in normal segments was 1.60 (1.29, 1.94) ml/min/g. The distribution of these two sets of values is shown in a boxplot in Fig. 5. The NMSE (standard deviation) between the predicted and measured AIFs was 0.002 (7.45).

Table 2

The median value (25th percentile, 75th percentile) of kinetic parameters estimates for the patient data, using the Combined PINN model.

Parameter	Median	(25th percentile, 75th percentile)
F_b (ml/min/g)	1.40	(1.12, 1.83)
v_b (dimensionless)	0.06	(0.05, 0.08)
v_e (dimensionless)	0.12	(0.10, 0.16)
PS (ml/min/g)	0.45	(0.30, 0.59)

Fig. 3

The estimated MBF maps for a patient (number 8) with no significant CAD (top) and patient (number 5) with CAD (bottom). The maps are shown under the corresponding MR images. The contrast of the MR images has been stretched to try to more clearly visualise the ischaemic regions.

Fig. 4

The training curves which show both the evolution of the loss terms (top two rows) and the (mean) estimated parameter values (bottom two rows) over the course of the training.

Fig. 5

The distribution of inferred MBF values in AHA segments with ischaemia (left) and without ischaemia (right).

Discussion

The aim of this study was to demonstrate the feasibility of inferring tracer-kinetic parameters from DCE-MRI data using physics-informed neural networks. PINNs are used to solve supervised learning tasks while respecting the underlying physics of the problem, usually represented in the form of differential equations. This allows the training of universal function approximators to solve differential equations that are highly data-efficient due to the enforcement of the physical laws. That is, PINNs do not need large databases of paired examples and labels, as are used in purely data-driven learning schemes. This paradigm can be extended beyond solving differential equations to inference problems which considers the kinetic parameters to be inferred as trainable parameters of the PINN.

Since first introduced by Raissi et al. (2019), PINNs have been employed for both solving differential equations and estimating parameters across a wide range of physical problems (Kadeethum, Jørgensen, Nick, 2020, Tartakovsky, Marrero, Perdikaris, Tartakovsky, Barajas-Solano, 2020), including cardiac electrophysiological modelling (Sahli Costabal et al., 2020) and blood flow modelling (Kissas et al., 2020). In this study, PINNs of different forms were developed and their performance was tested and compared to the traditional non-linear least squares fitting method in a simulated environment. In these synthetic experiments, the Combined PINN approach was shown to outperform the NLLS solutions based on the two quantitative metrics used: NMSE and SSIM. However, there are patch-like artifacts present in the PINN maps, showing the coupling between the kinetic parameter estimates and a limitation of the PINN approach. The DRO, as previously described by Debus et al. (2019), was chosen to minimise bias, but it is possible that a different DRO could lead to different results. The methods were further applied to patient data and yielded quantitative values close to the expected ranges. The MBF values are possibly lower than typically reported due to the high prevalence of ischaemia in these patients. The MBF values were also shown to discriminate well between normal and ischaemic tissue. This demonstrates the feasibility of the PINN-based kinetic parameter estimation in vivo, however, further work is required to properly validate the method, particularly in comparison to the established NLLS approach. To the best of our knowledge, this study presents the first use of PINNs for the inference of kinetic parameters in myocardial perfusion MR and DCE-MRI in general. Despite the fact that this is an initial proof-of-concept study, the promising results mean that it warrants further study which could lead to the adoption of PINNs as the method of choice for kinetic parameter estimation. Particularly since the adoption of quantitative MRI in clinical practice has been slow due to concerns about the accuracy and standardisation of methods.

Based on the DRO study, the Combined PINN performed the best of the presented methods with low NMSE and high SSIM as compared to the ground-truth kinetic parameters. This also shows the benefits of the flexibility of the PINN approach. It is seen that the Reduced PINN improves the results of the estimation as compared to the 2CXM PINN while the other three parameters suffer. However, a combination of the two models, in the Combined PINN, gives better results for all parameters. The different formulations of the models give different residuals that are optimised and these residuals have different weightings on different parameters. The Reduced model puts a heavier weighting on so it performs better for this parameter but not as well for the others. The 2CXM model and the Combined model more equally weight the parameters, and thus perform better on average across all parameters. is consistently the most difficult parameter to estimated and has high NMSE and low SSIM metrics. However, since the same trend is found with all methods it may be to do with the data sampling rather than the fitting. This is a well-known limitation of stress myocardial perfusion MR acquisitions that may not be long enough to measure the complete wash-out of the contrast agent and thus the full information regarding the kinetics in the extracellular-extravascular space may not be present in the data. Typical stress myocardial perfusion MR studies only report MBF values, however, as shown here, it is feasible to also assess the other kinetic parameters. It is possible that these parameters will give more in-depth insight to the microvasculature, a topic which is gaining increased attention in the literature (Rahman, Ryan, Lumley, Modi, McConkey, Ellis, Scannell, Clapp, Marber, Webb, Chiribiri, Perera, 2019, Rahman, Scannell, Demir, Ryan, McConkey, Ellis, Masci, Perera, Chiribiri, 2020).

While PINNs are based on tracer-kinetic models, they offer all of the flexibility deep learning and leverage the recent advances in the field. This includes the flexible optimisation algorithms and loss functions that may be combined and changed at will, and may be designed with intuition based on the governing physical law system. For example, in this study it was noted that the residual loss terms may be calculated for any time point , as the physical equations should hold regardless of the time input. So, while the observed measurements are fit at a limited number of points with a set temporal resolution (limited by the image acquisition process), the residuals can be enforced at any time point. The flexibility of the loss functions that can be optimised in deep learning frameworks was also shown and was used in this work to enforce non-negative concentration values. For this a continuous representation of the AIF is required and this is thus estimated by the neural network. The predicted AIFs well match the measured AIFs at the acquired time points, as shown by the low NMSE, with an example plotted in Fig. C1 and examples of voxel-wise fitted myocardial tissue curves in Fig. C2. Another property of conventional tracer-kinetic model fitting is that parameters are estimated pixelwise. That is, the estimates in each pixel are computed independently and spatial dependencies are rarely considered. However, these spatial correlations are exploited more easily in the PINN framework as a single neural network predicts parameters for all pixels simultaneously. Additionally, the spatial coordinates could be included as an input to the PINN to regularise the solution. The benefits of the joint optimisation of all kinetic parameters is clearly evident in Fig. 2 as the PINN approaches show much fewer outliers than the pixel-by-pixel fitting of the NLLS method.

Fig. C1

A comparison between a neural network predicted continuous time AIF and the discrete patient MRI measurements averaged from the LV, in terms of normalised signal intensity (SI) values.

Fig. C2

A comparison between two example neural network predicted myocardium tissue curves (with the decomposition into plasma and interstitial compartments) and the equivalent MRI measurements for a single myocardial voxel, in terms of normalised SI values, with the estimated kinetic parameters for each voxel also given.

While, in this work, different loss functions were tested and different inputs to the PINN were tried, there is also a huge potential to exploit the flexible framework further. For example, it would be possible to change the input from the image-derived concentration curves to the acquired k-t space and to learn both the image reconstruction and kinetic parameter estimation, similar to Dikaios (2020). However, this is also a potential limitation as there are a wide range of design choices and hyperparameters to be set, which have the potential to influence the results obtained. Recent work has focused on analysing the optimisation of PINNs (Wang et al., 2022) and the automatic choice of weighting parameters for the loss function (de Wolff et al., 2021). This work will likely serve to improve the accuracy and reproducibility of parameters estimated using PINNs.

It has been demonstrated previously that it is possible to train deep neural networks to directly estimate kinetic parameters (Scannell et al., 2019). However, this used a training set with corresponding labelled kinetic parameters. The difficulty of this approach is that there is not a ground-truth available for the training labels and these were derived from a separate inference scheme. The benefit of PINNs over such purely data-driven approaches is that it does not require labelled training data.

An important limitation of the proposed approach to consider for the inference of kinetic parameters is the computational cost as each individual slice or patient requires the training of a neural network. The proposed PINN method requires approximately 1 h per imaging slice utilizing an NVIDIA GeForce GTX 980M series GPU, as opposed to 3 min for the NLLS fitting. A possible solution to this could be the use of transfer learning. Using the assumption that most patients have similar kinetic parameters and concentration curves, a baseline PINN could be trained sequentially from a number of different patients and future patients could be processed by fine-tuning these weights rather than training from scratch. Initialising the weights from previously trained examples may also reduce spurious updates and increase robustness.

Conclusion

In this work, we have shown that it is feasible to perform myocardial perfusion quantification with physics-informed neural networks. Though the inference of kinetic parameters using PINNs is still slow and future work is required, this framework provides a high degree of flexibility and the initial results are promising. The methods for training neural networks are also still evolving at a high pace, making PINNs a good basis for future research. This may allow physics-informed machine learning to become an established method for tracer-kinetic modelling and quantitative MRI in general.

CRediT authorship contribution statement

Rudolf L.M. van Herten: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing – original draft, Writing – review & editing. Amedeo Chiribiri: Writing – review & editing, Supervision, Data curation, Project administration, Funding acquisition. Marcel Breeuwer: Conceptualization, Writing – review & editing, Supervision, Project administration, Funding acquisition. Mitko Veta: Conceptualization, Methodology, Writing – original draft, Writing – review & editing, Supervision, Project administration. Cian M. Scannell: Conceptualization, Methodology, Writing – original draft, Writing – review & editing, Supervision, Project administration.

Declaration of Competing Interest

Marcel Breeuwer is an employee of Philips Healthcare. All other authors declare that they have no known competing financial interests or personal relationships.