A Multi-Layer Perceptron Network for Perfusion Parameter Estimation in DCE-MRI Studies of the Healthy Kidney

Background: Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is an imaging technique which helps in visualizing and quantifying perfusion—one of the most important indicators of an organ’s state. This paper focuses on perfusion and filtration in the kidney, whose performance directly influences versatile functions of the body. In clinical practice, kidney function is assessed by measuring glomerular filtration rate (GFR). Estimating GFR based on DCE-MRI data requires the application of an organ-specific pharmacokinetic (PK) model. However, determination of the model parameters, and thus the characterization of GFR, is sensitive to determination of the arterial input function (AIF) and the initial choice of parameter values. Methods: This paper proposes a multi-layer perceptron network for PK model parameter determination, in order to overcome the limitations of the traditional model’s optimization techniques based on non-linear least-squares curve-fitting. As a reference method, we applied the trust-region reflective algorithm to numerically optimize the model. The effectiveness of the proposed approach was tested for 20 data sets, collected for 10 healthy volunteers whose image-derived GFR scores were compared with ground-truth blood test values. Results: The achieved mean difference between the image-derived and ground-truth GFR values was 2.35 mL/min/1.73 m2, which is comparable to the result obtained for the reference estimation method (−5.80 mL/min/1.73 m2). Conclusions: Neural networks are a feasible alternative to the least-squares curve-fitting algorithm, ensuring agreement with ground-truth measurements at a comparable level. The advantages of using a neural network are twofold. Firstly, it can estimate a GFR value without the need to determine the AIF for each individual patient. Secondly, a reliable estimate can be obtained, without the need to manually set up either the initial parameter values or the constraints thereof.


Introduction
Monitoring the state of kidney functioning has become an increasingly important task in modern medical diagnostics. Failures in renal operation may impair physiological homeostasis, leading to the improper management of electrolytes, acid-base balance perturbation, or the deregulation of arterial blood pressure. Kidneys are responsible for blood filtration, and the removal of water-soluble waste products of metabolism and surplus glucose and other organic substances. Kidney diseases include, among others, acute kidney injury [1], chronic kidney disease (CKD) [2] and various cancers (e.g., renal cell carcinoma). Treatment depends on the pathological condition, and may require life-long dialysis, or kidney removal or transplantation. In any case, precise knowledge about kidney performance is needed for the proper diagnosis, treatment and follow-up prognosis. Especially in the case of CKD, which may be caused by longstanding hypertension or diabetes mellitus, it is important to continuously and precisely monitor renal function, since early detection of the disease allows the prevention of its development to the end stage.
Tissue perfusion is a viable means of preserving cell nutrition, humoral communication and the elimination of waste products over a lifetime [3]. Thus information about local blood flow and the blood filtration taking part in the glomeruli plays a fundamental role in the diagnosis and follow-up of the aforementioned diseases. Routinely, renal diagnosis is accomplished through the creatinine clearance procedure, or by using the more specific iohexol clearance test. These tests lead to determination of the glomerular filtration rate (GFR)-the main indicator of kidney performance. However, neither procedure enables differentiation between left and right kidney function [4].
Dynamic contrast-enhanced (DCE) magnetic resonance imaging (MRI) is an alternative technique available for the assessment of kidney functioning [5][6][7]. The procedure requires the administration of a contrast agent (CA) intravenously, given as a bolus injection, whose passage through the abdominal arterial system and the capillary bed and tubular systems of the kidneys is measured. In contrast to conventional methods, DCE-MRI enables the estimation of GFR for a single kidney, ultimately at the voxel-level, and simultaneously it provides anatomical characteristics of the renal parenchyma.
The estimation of renal filtration based on the DCE-MR images involves the pharmacokinetic (PK) modeling of the contrast agent, arriving from the abdominal aorta into a kidney artery, the renal parenchyma and the afferent arterioles of the renal corpuscles, and then passing from the glomerular capillaries through the filtration barrier into the Bowman's space and the tubular system, ending up in the renal pelvis, the ureter and then the bladder as a constituent of the urine. There have been numerous PK models proposed in the literature [8], which vary in (1) the number of tissue compartments, (2) the level of incorporating compartment-specific impulse residue functions, (3) the inclusion of the effect of tracer agent leakage from capillaries in cancerous cells, and (4) the way in which the delay and the broadening of peak CA concentration in the renal parenchyma, relative to the abdominal aorta, is modeled. In the case of the normal-state kidney, all recent models exploit a two-compartment approach [9]. In this setting, at any point in time, the contrast agent concentration in the tissue is the sum of concentrations in the intravascular (IV) and extracellular extravascular (EEV) compartments. In the present study, we used the two-compartment filtration model (2CFM) proposed by Tofts et al. [10].
One approach to finding the PK model's parameters (defined below) is to fit the model-based signal intensity time course to the image data using a non-linear least-squares (NLLS) method. Depending on the model, various explicit parameters are used to characterize tracer agent kinetics. The so-called primary or independent parameters include plasma volume, plasma flow, interstitial volume and permeability product [9]. Based on these attributes, one may derive dependent parameters, i.e., the extraction fraction and the volume transfer constant. The latter, denoted as K trans , characterizes regional renal filtration, as it measures the ratio of plasma flow that is transferred from the glomerular capillary bed into the Bowman´s space. K trans multiplied by the kidney volume leads to an estimate of the glomerular filtration rate, and is therefore a central descriptor in quantifying kidney function.
In this study, we propose to use a multi-layer perceptron (MLP) network for PK model parameter estimation, instead of the non-linear least-squares curve-fitting method. The rationale behind this concept arises from the fact that the results of the iterative methods may be located in a local and-in general-suboptimal solution. Thus, the outcome depends largely on the initial choice of parameter values. Moreover, one usually has to establish constraints on the expected parameter values, so that the result falls into a physiologically feasible range. Bounds that are too narrow, however, may prevent the fitting algorithm from finding an optimum. Contrary to NLLS, MLP does not require that one provides the initial parameter values; moreover, it ensures the optimal solution when properly trained.
Although MLP networks could be regarded as a legacy method, they are still in use, and prove efficient in complex non-linear regression problems [11,12]. Moreover, according to the universal approximation theorem, an artificial neural network with at least one hidden layer of nonpolynomial activation (such as sigmoid or rectified linear) can model any measurable function that maps one finite-dimensional space to another, provided the hidden layers contain enough units [13]. Therefore, MLP seems to be particularly usefully in its application to the specific problem of PK model parameter estimation in renal DCE-MRI examinations.
Also, the reliability of the fitted PK model parameters is conditioned by-among other factors-the determination of the arterial input function (AIF), i.e., the tracer concentration's time course in the feeding artery. The standard way of annotating AIF consists of estimating it for each patient. However, other approaches are accepted if patient-specific AIF measurements are not possible. Such problems can arise due to data acquisition constraints or the lack of a suitable artery within the imaging field of view. In the alternative strategy, one either utilizes a population-based mean AIF [14] (calculated for a cohort of studies, where it is feasible to do so) or defines its parameterized functional form [15]. In this study, we implicitly follow the population-based approach. The neural network is trained using a pool of AIFs automatically derived for a subset of available DCE images. Then, the GFR for an individual test case is estimated without the need to determine the case-specific AIF.
To sum up, the main contribution of this paper is the establishment of the architecture of an MLP network for estimating parameters of the 2CFM, based solely on image intensity time courses in the kidney region. Additionally, a side-contribution is the introduction of the algorithm for automatic determination of the AIF in DCE-MRI examinations. The latter achievement was needed in order to construct the data set for MLP training, but the algorithm can be equally used in standard approaches to PK modeling.

Subjects
A total of 20 DCE-MRI examinations were available for experiments [16]. The data sets were collected from 10 healthy volunteers. Each subject was imaged twice, seven days apart (further, these examinations will be referred to as Session 1 and 2). The acquisition sequence used the standard 3D fast low-angle shot (FLASH) spoiled gradient recalled echo technique with the following parameters: TR = 2.36 ms, TE = 0.8 ms, FA = 20 • , in-plane resolution = 2.2 × 2.2 mm 2 , slice thickness = 3 mm, acquisition matrix = 192 × 192, number of slices = 30. Prior to image acquisition, patients were administered 0.025 mmol/kg of GdDOTA at 3 mL/s flow rate. The contrast agent was injected intravenously. Then, 74 volumetric scans were gathered in 2.3 s time intervals.
Every image series was registered in the time domain using the b-splines method, as implemented in the Insight Toolkit (ITK) library [17]. Specifically, we used Slicer 3D software to accomplish the task. For every study, the same reference frame was used as a fixed volume that every other (moving) volume was matched to. B-splines registration was performed fully automatically, i.e., no fiducial points were marked over tissues of interest. Furthermore, the procedure was launched in a multi-stage configuration. In each stage, various settings of grid size and subsampling rates were used. For a detailed interpretation of these parameters, the reader is referred to the ITK documentation. In short, they allow the registration of images at different scales-starting from coarse matching and then refining the outcome.
A common practice in DCE series analysis is to interpolate the measured signal so as to achieve a higher temporal resolution and capture more precisely the critical phases of the perfusion process, such as, e.g., bolus arrival time or the peak of the signal enhancement [10]. As such, we applied cubic spline interpolation to the observed intensity time courses, and then resampled the interpolated signal curves at 0.4 milliseconds intervals.
As stated in the introduction, analysis of the DCE images leads to estimation of the GFR. Thus, in order to validate the obtained estimates against ground-truth values, volunteers underwent iohexol clearance tests. The measurement was carried out by administrating a dose of 5 mL of iohexol (300 mg I/mL; Omnipaque 300, GE Healthcare), and then acquiring a venous blood sample after 4 h.
No specific diet was imposed on the participating subjects, but they were asked to abstain from alcohol and high-protein meals, avoid exhausting physical exertion, be normally hydrated at least 2 days before examination, and have no caffeine on the examination day. Apart from these restrictions, participants were instructed to retain their ordinary habits concerning the type and amount of consumed food so as to assure stable examination conditions. All subjects gave their informed consent for inclusion before they participated in the study. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Regional Committees for Medical Research Ethics Western Norway (REC West 2012/1869).

Two-Compartment Filtration Model of Kidney Perfusion
In the PK model considered in this study [10], it is assumed that renal tissue is composed of two compartments-an intravascular (IV) one and an extracellular extravascular (EEV) one ( Figure 1). The model can be applied either to the whole kidney parenchyma or the kidney cortex only.
The assumption which could be violated here is that CA does not flow out from the tubular system within the fitting period.

Arterial Input Function
Although in the proposed approach, the trained neural network enables the determination of GFR without annotation of a patient-specific arterial input function, information about the time course of the contrast agent's concentration in a tissue-feeding artery is required at the training stage. As further described, we use a subset of image-derived AIFs to construct a training data set. Therefore, in this section, we present our custom method for automatic AIF determination.
Most frequently, AIF is determined based on signal averaged over a manually delineated region of interest near to the perfused organ [18,19]. Specifically, in the case of renography it is recommended to locate the AIF region-of-interest (ROI) near the place where the aorta bifurcates into common iliac arteries. This location is in relative proximity to the kidney, and is simultaneously sufficiently far from the upper aorta section, where the signal coming from the blood that enters the imaging slab at high velocity is high, because of the inflow artifact and not because of the contrast agent [19]. Moreover, Cutajar et al. in [18] confirmed the significant dependence of renal perfusion and the filtration coefficients on the size of the chosen AIF ROI. In [20], AIF was determined from voxels chosen semi-automatically on the maximum signal-enhancement maps as the brightest locations in the abdominal aorta, directly below the inlets of the renal arteries. Note that the manual or semi-automatic approaches suffer from poor reproducibility and subjectivism, and they can be time-consuming and can require the availability of well-trained staff. In addition, it may be necessary in such cases to obtain patient-specific hemodynamic parameters in order to adjust the DCE-MRI protocol to the arterial pulse wave, and thus reduce the inflow artifact. Therefore, there were several The contrast agent's concentration in the renal tissue C tissue (t) is governed by the equation where C art p denotes the arterial input function, v p -plasma volume fraction, and C kid p -CA concentration in the plasma, given by the convolution The first term in (1) represents the EEV compartment, whereas Equation (2) models the delay and dispersion of the arterial input upon its arrival to the IV compartment by using the concept of the vascular impulse response function (VIRF), denoted as g(t). In our study we used the delayed exponential VIRF, defined as The variables T g -the dispersion time constant-and ∆-the delay interval-together with the volume fraction v p and transfer constant K trans form the complete set of the Tofts model's parameters.

Arterial Input Function
Although in the proposed approach, the trained neural network enables the determination of GFR without annotation of a patient-specific arterial input function, information about the time course of the contrast agent's concentration in a tissue-feeding artery is required at the training stage. As further described, we use a subset of image-derived AIFs to construct a training data set. Therefore, in this section, we present our custom method for automatic AIF determination.
Most frequently, AIF is determined based on signal averaged over a manually delineated region of interest near to the perfused organ [18,19]. Specifically, in the case of renography it is recommended to locate the AIF region-of-interest (ROI) near the place where the aorta bifurcates into common iliac arteries. This location is in relative proximity to the kidney, and is simultaneously sufficiently far from the upper aorta section, where the signal coming from the blood that enters the imaging slab at high velocity is high, because of the inflow artifact and not because of the contrast agent [19]. Moreover, Cutajar et al. in [18] confirmed the significant dependence of renal perfusion and the filtration coefficients on the size of the chosen AIF ROI. In [20], AIF was determined from voxels chosen semi-automatically on the maximum signal-enhancement maps as the brightest locations in the abdominal aorta, directly below the inlets of the renal arteries. Note that the manual or semi-automatic approaches suffer from poor reproducibility and subjectivism, and they can be time-consuming and can require the availability of well-trained staff. In addition, it may be necessary in such cases to obtain patient-specific hemodynamic parameters in order to adjust the DCE-MRI protocol to the arterial pulse wave, and thus reduce the inflow artifact. Therefore, there were several attempts made to automate determination of the AIF. A common approach employs clustering techniques [21]. On the other hand, in [22], the authors used Kendall's coefficient of concordance. These attempts, however, were either applied in preclinical studies, or to different organs such as the head or neck, but not to the kidney.
We postulate a fully automatic procedure for determining a patient-specific AIF based on the clustering of abdominal aorta voxels' time-series. Our approach consists of three stages, as illustrated in Figure 2. The algorithm starts with delineation of the abdominal aorta region. This task, further decomposed into four steps, is accomplished through the vessel segmentation of one time frame of the DCE-MR image series. The appropriate time frame is selected for maximum signal enhancement in the bottom-central region of the image, near the aorta's bifurcation into common iliac arteries. The peak intensity value in this region corresponds to the time point when the CA fills the whole aorta lumen, thus ensuring the highest contrast between the background tissue and the aorta at its entire length.
Next, the image is blurred by the low-pass Gaussian filter (σ = 3 voxels) to suppress fine structures that could be misinterpreted as vessels. Such a smoothed image is then submitted to the vessel enhancement routine. In this step, we computed a multi-scale vesselness function in the form proposed by Sato [23], using 5 scales of radius standard deviation ranging from 0.5 to 2.5 mm. As a result, the filtered image contains only the elongated bright structures. Apart from the aorta, it may contain smaller vessels, like renal arteries and some other elongated structures. In order to remove them, the flood fill operation is executed, with the lower threshold set to an empirically found value equal to 12.5% of the maximum image intensity. The upper bound is set to maximum, as we found the aorta to be the brightest structure visible in the time frames selected for segmentation. Similarly, the seed point of the flood fill operation was determined as the brightest point in the central part of the image along the coronal direction.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 22 contain smaller vessels, like renal arteries and some other elongated structures. In order to remove them, the flood fill operation is executed, with the lower threshold set to an empirically found value equal to 12.5% of the maximum image intensity. The upper bound is set to maximum, as we found the aorta to be the brightest structure visible in the time frames selected for segmentation. Similarly, the seed point of the flood fill operation was determined as the brightest point in the central part of the image along the coronal direction. In the second stage, abdominal aorta voxels are grouped into clusters using their signal intensity time courses as the feature vectors. For this purpose, we employed a k-means clustering algorithm with k = 5. We found this value to be a reasonable tradeoff between the need to reflect apparent diversity in the aorta signal dynamics and the need to keep the cluster size sufficiently large, so that each one represents a significant portion of aorta voxels. As an effect of clustering, points whose temporal signal characteristics are similar are joined. Eventually, we select one of the clusters as the AIF ROI. The selection is based on the minimum mean signal value calculated in the pre-bolus phase. We assume that the best AIF candidate region has the lowest signal in this initial phase of bolus passage, as it is free from the inflow-enhancement artifact.

PK Model Fitting
There are several optimizers utilized in the NLLS curve-fitting problems. Downhill Simplex, Levenberg-Marquardt or non-linear conjugate gradients belong to the most frequently cited algorithms [24]. In the DCE-MRI context, a popular solver is based on the gradient expansion method implemented in the IDL Software [7]. Often, in order to ensure that the resulting parameter set falls into a range of reasonable physiological values, it is also necessary to restrain the optimization procedure to a specific solution subspace. In such a case, the Trust Region Reflective technique can be considered as the algorithm of choice as it enables optimization with constraints. All the abovementioned methods are iterative, and require initiation with a starting point representing a feasible solution. As such, there is a risk of getting stuck in a local minimum of the objective function, usually defined as the mean squared error of the model residuals.
Therefore, we propose to solve the curve-fitting problem using an artificial neural network (ANN); specifically, a fully-connected multi-layer perceptron structure. ANNs have been In the second stage, abdominal aorta voxels are grouped into clusters using their signal intensity time courses as the feature vectors. For this purpose, we employed a k-means clustering algorithm with k = 5. We found this value to be a reasonable tradeoff between the need to reflect apparent diversity in the aorta signal dynamics and the need to keep the cluster size sufficiently large, so that each one represents a significant portion of aorta voxels. As an effect of clustering, points whose temporal signal characteristics are similar are joined. Eventually, we select one of the clusters as the AIF ROI. The selection is based on the minimum mean signal value calculated in the pre-bolus phase. We assume that the best AIF candidate region has the lowest signal in this initial phase of bolus passage, as it is free from the inflow-enhancement artifact.

PK Model Fitting
There are several optimizers utilized in the NLLS curve-fitting problems. Downhill Simplex, Levenberg-Marquardt or non-linear conjugate gradients belong to the most frequently cited algorithms [24]. In the DCE-MRI context, a popular solver is based on the gradient expansion method implemented in the IDL Software [7]. Often, in order to ensure that the resulting parameter set falls into a range of reasonable physiological values, it is also necessary to restrain the optimization procedure to a specific solution subspace. In such a case, the Trust Region Reflective technique can be considered as the algorithm of choice as it enables optimization with constraints. All the above-mentioned methods are iterative, and require initiation with a starting point representing a feasible solution. As such, there is a risk of getting stuck in a local minimum of the objective function, usually defined as the mean squared error of the model residuals.
Therefore, we propose to solve the curve-fitting problem using an artificial neural network (ANN); specifically, a fully-connected multi-layer perceptron structure. ANNs have been successfully employed in a variety of optimization scenarios, including the fitting of complex, multi-parametric relationships to the observed noisy data [25]. The neural network training process is configured to explore a wide range of feasible parameter values, potentially covering the whole solution space. For every parameter candidate set, there is a corresponding CA concentration time course generated based on the PK model equation. Thus, the network incorporates a priori knowledge about the underlying mechanism of signal generation. Then, the network's response to an observed, patient-specific CA time-course is estimated in the generalized context of tracer kinetics, making the fitting process more robust to noise, measurement errors and local minima.
The first step in the design of an MLP network is the choice of the number of hidden layers. As formulated by the already recalled universal approximation theorem, one sufficiently large hidden layer should enable the accurate approximation of the modeled function. In cases of complex problems, however, two layers with lower numbers of units may be easier to train, and could better acquire a general pattern that links input data vectors with estimated function parameters. We have observed such behavior in the current study. One-hidden-layer structures turned out to be more prone to overfitting, and their training process was less liable to converge. Therefore, in the rest of the paper we concentrate on the description of the three-layer architecture, i.e., consisting of two hidden layers and one output layer. However, in the Results section we provide example outcomes of training two-layer structures to support the above observations. On the other hand, increasing the network depth to three hidden layers did not improve training efficacy, while making the task of optimizing the network architecture more complicated. Hence, the results of these computations are omitted in the experiments description.
The overall concept of a three-layer perceptron network operating as an estimator of the Tofts model parameters is shown in Figure 3. The network design is represented by weight vectors w = [w 1 , w 2 ], and v for two hidden layers and the output layer neurons respectively. The symbols f (·) and h(·) denote the activation functions of perceptrons in a given layer. In our implementation we used the same activation for both hidden layers, which was the rectified linear unit function, i.e., h(x) = max(x, 0). For the output layer we used linear activation. The detailed structure of implemented MLP is presented in the Results section.
In principle, the network inputs accept a vector C tissue (t) of CA concentrations in the subsequent time frames of a given DCE series. Transformation from the image signal intensity S(t) to concentration C tissue (t) is given in the Appendix A. In the following, the time variable t will be skipped where possible, for the sake of clarity.
There are two modes of network operation. In the recall mode, a trained network recognizes the pattern in the observed data and predicts the possible values of the PK model's parameters. Formally, the network realizes the transformationθ = f (v, h(w, C tissue )), which approximates the unknown mappingθ = F (C tissue ) from the observation space to the parameter space. It is assumed that this mapping exists and is unique. In the training mode, the network is fed by a large and versatile collection of input vectors, and simultaneously it is presented with an expected parameter set θ for each individual input example. The difference between these true parameter values and the output drives the process of network weights adjustment. Specifically, we applied the mean squared error (MSE) metric as the network loss function. Apart from MSE, during the network weights adjustment, the mean absolute percentage error (MAPE) was also calculated for both the training and validation samples. The definitions of the MSE and MAPE metrics are given in Appendix A.
The weights assigned to random values are updated after each iteration using the error backpropagation algorithm [26]. The loss function (dependent on inter-layer weights) was optimized using the stochastic gradient descent algorithm [27]. The network was implemented in a Python script using Keras library with the Tensorflow backend [28].
The overall concept of a three-layer perceptron network operating as an estimator of the Tofts model parameters is shown in Figure 3. The network design is represented by weight vectors w = [w1, w2], and v for two hidden layers and the output layer neurons respectively. The symbols f(•) and h(•) denote the activation functions of perceptrons in a given layer. In our implementation we used the same activation for both hidden layers, which was the rectified linear unit function, i.e., h(x) = max(x, 0). For the output layer we used linear activation. The detailed structure of implemented MLP is presented in the Results section.  All layers are fully connected, as in the conventional multi-layer perceptron (MLP) structures, except dropout mechanism was switched on, in which case random connections between hidden layers were removed during training to prevent network overfitting.
In principle, the network inputs accept a vector Ctissue(t) of CA concentrations in the subsequent time frames of a given DCE series. Transformation from the image signal intensity S(t) to concentration Ctissue(t) is given in the Appendix A. In the following, the time variable t will be skipped where possible, for the sake of clarity.
There are two modes of network operation. In the recall mode, a trained network recognizes the pattern in the observed data and predicts the possible values of the PK model's parameters. Formally, the network realizes the transformation = , , , which approximates the unknown mapping = from the observation space to the parameter space. It is assumed that this mapping exists and is unique. In the training mode, the network is fed by a large and versatile collection of input vectors, and simultaneously it is presented with an expected parameter set for each individual input example. The difference between these true parameter values and the output drives the process of network weights adjustment. Specifically, we applied the mean squared error (MSE) metric as the network loss function. Apart from MSE, during the network weights adjustment, the mean absolute percentage error (MAPE) was also calculated for both the training and validation samples. The definitions of the MSE and MAPE metrics are given in Appendix A.
The weights assigned to random values are updated after each iteration using the error backpropagation algorithm [26]. The loss function (dependent on inter-layer weights) was optimized using the stochastic gradient descent algorithm [27]. The network was implemented in a Python script using Keras library with the Tensorflow backend [28].
The crucial step in the above-described procedure is in the preparation of the training set. On the one hand it must be sufficiently large and varied for the network to gain the generalization capabilities. On the other hand, the input signal intensity vectors must be accompanied by the true parameter values, which allow the observation of a given Ctissue in response to a specific arterial input function. Moreover, since the available data set size is relatively small, we repeat this procedure using the leave-one-out approach. Having a data set of 20 examinations, we invoked calculations 20 times, Figure 3. Overview of a pharmacokinetic (PK) model parameter fitting using an artificial neural network operating in the training (a) and recall (b) modes. The three-layer perceptron structure is shown only conceptually-the number of hidden units was adjusted during the network design process. The number of output units was always set to 4, i.e., the number of PK model parameters. All layers are fully connected, as in the conventional multi-layer perceptron (MLP) structures, except dropout mechanism was switched on, in which case random connections between hidden layers were removed during training to prevent network overfitting. The crucial step in the above-described procedure is in the preparation of the training set. On the one hand it must be sufficiently large and varied for the network to gain the generalization capabilities. On the other hand, the input signal intensity vectors must be accompanied by the true parameter values, which allow the observation of a given C tissue in response to a specific arterial input function. Moreover, since the available data set size is relatively small, we repeat this procedure using the leave-one-out approach. Having a data set of 20 examinations, we invoked calculations 20 times, each time excluding another subject AIF from the training stage. Hence, a training set in a given repetition was established via the 3-step procedure described below.
Step 1. Selecting a subset of 19 realistic arterial input functions.
Step 2. Establishing ranges of model parameter values. Table 1 presents the assumed parameter limits along with their corresponding probability distributions. These limits were fixed with respect to the published values for normal subjects [7,10,29]. Training input samples were generated using the model equation and parameters sampled from these distributions. According to the Tofts model [10], there are four such parameters: K trans , v p , T g and ∆. The latter two are the parameters of the vascular impulse response function, whose sum defines the so-called mean residence time (MRT). Since [7,10] report only the value of MRT, we sampled T g and MRT, whereas ∆ was calculated as Step 3. For a given AIF, we sample the model parameters θ and calculate C tissue according to Equation (1). Sampling is repeated 10 4 times for each AIF. As a result, there were 19 × 10 4 < C(t), θ > pairs available in each leave-one-out repetition. The calculated C tissue curves were probed at time steps adjusted to match the temporal resolution of the network input vector. The simulated tracer concentration time curves were optionally corrupted by the Gaussian noise to increase their variability and make them more realistic. Data corruption was performed by adding to each simulated CA concentration a value drawn from the normal distribution, with mean = 0, and standard deviation was adjusted to a current time step, i.e., σ noise = s·[C tissue (t) − min(C tissue (t))]. The scale factor s was experimentally set to 0.025 to achieve reasonable deflection from the modeled CA concentration time curves. Then, the data set was partitioned into training and test sets in the proportion 7:3. Moreover, 30% of the training data vectors were randomly set apart in each epoch for validation purposes.
The last issue which needs to be fixed is the establishing of a proper network architecture. During experiments, it occurred that only one hidden layer cannot accurately encapsulate the non-linear relationship between the PK model's parameters and the variety of CA concentration time-curves. The two-layer structure ensured a significant decrease in the loss function value, while the usage of more than two hidden layers did not offer any further improvement in this respect. Eventually, the number of perceptrons in each layer had to be established. We tested 12 configurations, listed in Table 2. Apart from the various cardinalities of neurons in hidden layers, we also considered inclusion of additional dropout layers between them to overcome the problem of overfitting [30]. For every configuration and each leave-one-out fold, the training was run for a fixed number of 100 epochs.

Results
We compared the effectiveness of the postulated ANN-based approach to optimizing a PK model's parameters in reference to the algebraic curve-fitting method, configured to utilize the Trust Region Reflective algorithm [31]. Apart from the direct comparison of particular model parameters, we also estimated single-and two-kidney GFRs. The image-derived estimates were then compared against reference values measured with iohexol clearance tests.
We also validated the designed algorithm for the determination of the arterial input function. In this experiment, we calculated the GFR scores using the Trust Region Reflective method for AIFs annotated manually, and using the proposed automated approach.
For every study we performed manual segmentation of the whole kidney, as well as cortex and medulla. This step was accomplished by two professionals experienced in medical image analysis, using our custom software annotation tool. The rate of agreement between the respective kidney regions delineated by both experts was estimated by the Jaccard coefficient, and was equal to 0.93 on average. For a given study, kidney segments were outlined only in a time frame corresponding to the maximum signal enhancement in the perfusion phase, which ensured the largest contrast between cortex, medulla and pelvis. Depending on the study, this maximal enhancement was observed in frames numbered from 12 to 18. Then, the segmentations were applied to all other time frames reflected in the parameter estimation procedure, so as to determine mean signal time courses in the perfused renal tissue, and the selected frames became fixed reference frames in the b-spline registration algorithm. Figure 4 visualizes the locations of the arterial input function ROIs, as identified by the proposed automatic method and manual selection. The presented examples correspond to the two cases, where the best and the worst model fits were obtained. As shown, in comparison to manual AIF annotations, automatically found ROIs occupy many more voxels that are free from inflow artifacts. As a result, fitted intensity time curves better encapsulate image signal changes. For example, in Figure 4c,e, the curve is visible in the region of maximum signal enhancement. In the case of automatically determined AIF, the curve approximates well the measured signal samples, whereas in the case of manual annotations, the model overestimates signal values. The reliability of the GFR measurements based on automatically determined AIFs is visualized in Figure 5b, and can be compared to the Bland-Altman plots shown in Figure 5a, obtained with manually annotated AIFs. Since the automatically determined AIFs seemed to improve the accuracy and precision of GFR assessment, this allowed us to use them in the generation of the training set for the neural network.

Assessment of Automatic AIF Determination Method
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 22

Results
We compared the effectiveness of the postulated ANN-based approach to optimizing a PK model's parameters in reference to the algebraic curve-fitting method, configured to utilize the Trust Region Reflective algorithm [31]. Apart from the direct comparison of particular model parameters, we also estimated single-and two-kidney GFRs. The image-derived estimates were then compared against reference values measured with iohexol clearance tests.
We also validated the designed algorithm for the determination of the arterial input function. In this experiment, we calculated the GFR scores using the Trust Region Reflective method for AIFs annotated manually, and using the proposed automated approach.
For every study we performed manual segmentation of the whole kidney, as well as cortex and medulla. This step was accomplished by two professionals experienced in medical image analysis, using our custom software annotation tool. The rate of agreement between the respective kidney regions delineated by both experts was estimated by the Jaccard coefficient, and was equal to 0.93 on average. For a given study, kidney segments were outlined only in a time frame corresponding to the maximum signal enhancement in the perfusion phase, which ensured the largest contrast between cortex, medulla and pelvis. Depending on the study, this maximal enhancement was observed in frames numbered from 12 to 18. Then, the segmentations were applied to all other time frames reflected in the parameter estimation procedure, so as to determine mean signal time courses in the perfused renal tissue, and the selected frames became fixed reference frames in the b-spline registration algorithm. Figure 4 visualizes the locations of the arterial input function ROIs, as identified by the proposed automatic method and manual selection. The presented examples correspond to the two cases, where the best and the worst model fits were obtained. As shown, in comparison to manual AIF annotations, automatically found ROIs occupy many more voxels that are free from inflow artifacts. As a result, fitted intensity time curves better encapsulate image signal changes. For example, in Figures 4c,e, the curve is visible in the region of maximum signal enhancement. In the case of automatically determined AIF, the curve approximates well the measured signal samples, whereas in the case of manual annotations, the model overestimates signal values. The reliability of the GFR measurements based on automatically determined AIFs is visualized in Figure 5b, and can be compared to the Bland-Altman plots shown in Figure 5a, obtained with manually annotated AIFs. Since the automatically determined AIFs seemed to improve the accuracy and precision of GFR assessment, this allowed us to use them in the generation of the training set for the neural network.     (a) (b) Figure 5. Bland-Altman plots of agreement for manually (a) and automatically (b) determined AIF ROIs. Measurements were evaluated against normality using Shapiro-Wilk test. The obtained pvalues were equal to 0.52 (iohexol-based glomerular filtration rate (GFRs)), 0.29 (image-derived GFRs with manual ROIs) and 0.12 (for image-derived GFRs with automatically determined ROIs), thus disallowing us to reject the null hypothesis of the data with a normal distribution.

Selecting Network Configuration
In the first phase of experiments, we aimed to select an architecture for the neural network suitable to solving the problem of fitting the 2CFM to the DCE-MRI data. We began with deciding upon the number of hidden layers and their activation units. Figure 6 presents the learning curves for three configurations of two-layer architectures, whereas the plots shown in Figure 7 correspond to the three-layer network configurations listed in Table 2. Although the selected learning curves were obtained only for one of the leave-one-out repetitions, they are representative for all the subjects in the study. It can be seen that in all cases, the loss function initially decreases rapidly, and converges to a value in the range 0.01-0.13 depending on the network architecture and the existence of noise in the training data. However, for configurations without the dropout layer, it is observable that the MSE value in a validation subset heavily fluctuates until the very end of the training process. It must be noted that for the single-hidden-layer configurations (Figure 6), the problem with the convergence of the loss function also persists in the largest tested network, with 50 hidden units, despite even inserting a dropout layer before the output nodes. This observation, combined with potential overfitting (which manifests in larger errors in a validation sample, with respect to a training set), substantiates the use of two hidden layers. In the latter case, the loss function stabilizes after around 40 iterations for 'dropout' configurations with higher numbers of neurons (architectures number 6, 8, 10 and 12). It is especially pronounced if the training data was corrupted by noise.

Selecting Network Configuration
In the first phase of experiments, we aimed to select an architecture for the neural network suitable to solving the problem of fitting the 2CFM to the DCE-MRI data. We began with deciding upon the number of hidden layers and their activation units. Figure 6 presents the learning curves for three configurations of two-layer architectures, whereas the plots shown in Figure 7 correspond to the three-layer network configurations listed in Table 2. Although the selected learning curves were obtained only for one of the leave-one-out repetitions, they are representative for all the subjects in the study. It can be seen that in all cases, the loss function initially decreases rapidly, and converges to a value in the range 0.01-0.13 depending on the network architecture and the existence of noise in the training data. However, for configurations without the dropout layer, it is observable that the MSE value in a validation subset heavily fluctuates until the very end of the training process. It must be noted that for the single-hidden-layer configurations (Figure 6), the problem with the convergence of the loss function also persists in the largest tested network, with 50 hidden units, despite even inserting a dropout layer before the output nodes. This observation, combined with potential overfitting (which manifests in larger errors in a validation sample, with respect to a training set), substantiates the use of two hidden layers. In the latter case, the loss function stabilizes after around 40 iterations for 'dropout' configurations with higher numbers of neurons (architectures number 6, 8, 10 and 12). It is especially pronounced if the training data was corrupted by noise.  Figure 6. Learning curves for various two-layer MLP configurations-training set scores in blue, validation set scores in orange. The plots in the first two and second two columns correspond to architectures without and with the dropout layer, respectively. The curves in columns 1 and 3 were obtained for the noise-free training data, whereas in columns 2 and 4 for the data corrupted with Gaussian noise.
Apparently, increased data variability due to the presence of noise, along with the inclusion of a dropout layer, ensures better loss function convergence and helps prevent the network from overfitting. As can be observed in Figure 7 (second column), the loss function scores on the validation set start to surpass the values obtained for the training set quite early in the learning process. This overfitting dynamic is further confirmed by the MAPE scores obtained after 100 learning epochs (see Figure 6. Learning curves for various two-layer MLP configurations-training set scores in blue, validation set scores in orange. The plots in the first two and second two columns correspond to architectures without and with the dropout layer, respectively. The curves in columns 1 and 3 were obtained for the noise-free training data, whereas in columns 2 and 4 for the data corrupted with Gaussian noise. Apparently, increased data variability due to the presence of noise, along with the inclusion of a dropout layer, ensures better loss function convergence and helps prevent the network from overfitting. As can be observed in Figure 7 (second column), the loss function scores on the validation set start to surpass the values obtained for the training set quite early in the learning process. This overfitting dynamic is further confirmed by the MAPE scores obtained after 100 learning epochs (see Figure 8). Without the dropout mechanism, the error remains higher for the validation set than for the training one, for almost every network configuration.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 22 Figure 8). Without the dropout mechanism, the error remains higher for the validation set than for the training one, for almost every network configuration. Learning curves for various three-layer MLP configurations-training set scores in blue, validation set scores in orange. The plots in the first two and two second two columns correspond to architectures without and with the dropout layer, respectively. The curves in columns 1 and 3 were obtained for the noise-free training data, whereas in columns 2 and 4 for the data corrupted with Gaussian noise.
While the above analysis settles upon the usage of the dropout layer and the noisy training vectors, the question regarding the sizes of the hidden layers must still be resolved. On one side,  Learning curves for various three-layer MLP configurations-training set scores in blue, validation set scores in orange. The plots in the first two and two second two columns correspond to architectures without and with the dropout layer, respectively. The curves in columns 1 and 3 were obtained for the noise-free training data, whereas in columns 2 and 4 for the data corrupted with Gaussian noise.
predicated outcomes do not conform to the actual data distribution. This tendency is not observed in the case of noisy data; however, changing the number of neurons in the first hidden layer from 40 to 50 did not bring significant improvement in the R 2 score (0.85 versus 0.86), and simultaneously, as shown in Figure 8d, it slightly increased the mean absolute percentage error for the validation set (from 13.9% to 14.3%). Therefore, for the rest of the experiments we chose configuration #10, with 40 and 20 perceptrons in the first and second hidden layers, respectively, and with a dropout layer in between, with a dropout rate = 0.2. This setting ensures the lowest mean absolute percentage error, relatively high R 2 scores for all estimated parameters (Figure 10), and that the network model's capacity is appropriate for the data set's complexity. While the above analysis settles upon the usage of the dropout layer and the noisy training vectors, the question regarding the sizes of the hidden layers must still be resolved. On one side, increasing the number of neurons leads to the lowering of the MSE and MAPE scores. On the other hand, however, it escalates the risk of overfitting. Figure 9 presents the distributions of the ∆ parameter estimates against the true values in the test set. Regression analysis shows good agreement between the predicted and the actual parameter estimates. The performed t-test and its resulting large p-values indicate that one cannot reject the hypothesis of equality between the group means. However, the overfitting effect manifests itself in the case of noise-free data, and it appears deeper for the more complex MLP architectures. The network tends to predict parameter values close to some discrete levels. Although the determination coefficients R 2 are higher in these cases, clearly the predicated outcomes do not conform to the actual data distribution. This tendency is not observed in the case of noisy data; however, changing the number of neurons in the first hidden layer from 40 to 50 did not bring significant improvement in the R 2 score (0.85 versus 0.86), and simultaneously, as shown in Figure 8d, it slightly increased the mean absolute percentage error for the validation set (from 13.9% to 14.3%).
Therefore, for the rest of the experiments we chose configuration #10, with 40 and 20 perceptrons in the first and second hidden layers, respectively, and with a dropout layer in between, with a dropout rate = 0.2. This setting ensures the lowest mean absolute percentage error, relatively high R 2 scores for all estimated parameters (Figure 10), and that the network model's capacity is appropriate for the data set's complexity. Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 22

Evaluation of MLP-Based Model Fitting Reliability
The second phase of experiments consisted of using the trained MLP network to estimate 2CFM model parameters, as well as the GFR values for the 20 DCE-MRI examinations available in this study. The statistics of the parameter estimates are summarized in Table 3, whereas the calculated GFR values are collected in Table 4. The reliability of GFR estimation was also evaluated against the ground-truth using the Bland-Altman plot, shown in Figure 11. When compared to Figure 5b, it emerges that the mean biases from the reference values are similar for both optimization methods-MLP and NLLS curve-fitting. However, the limits of agreement and confidence intervals are narrower for the neural network. Moreover, our ANNs ensure the better balance of the estimated values around the mean. The NLLS method, on the other hand, tends to overestimate the measured quantity, especially in the lower range of values. On the other hand, in the case of ANN, GFR was estimated without explicit determination of AIF for a given subject. It proves that the designed network was capable of extracting the general characteristics of the underlying perfusion process.

Evaluation of MLP-Based Model Fitting Reliability
The second phase of experiments consisted of using the trained MLP network to estimate 2CFM model parameters, as well as the GFR values for the 20 DCE-MRI examinations available in this study. The statistics of the parameter estimates are summarized in Table 3, whereas the calculated GFR values are collected in Table 4. The reliability of GFR estimation was also evaluated against the ground-truth using the Bland-Altman plot, shown in Figure 11. When compared to Figure 5b, it emerges that the mean biases from the reference values are similar for both optimization methods-MLP and NLLS curve-fitting. However, the limits of agreement κ and confidence intervals are narrower for the neural network. Moreover, our ANNs ensure the better balance of the estimated values around the mean. The NLLS method, on the other hand, tends to overestimate the measured quantity, especially in the lower range of values. On the other hand, in the case of ANN, GFR was estimated without explicit determination of AIF for a given subject. It proves that the designed network was capable of extracting the general characteristics of the underlying perfusion process.
From the analysis of Table 3, it emerges that MLP determined K trans and v p parameter values close to those of the reference method. The box plots presented in Figure 12a,b visually illustrate that the ranges of estimated K trans and v p values partially overlap, and their medians are similar. The result of the t-test gives evidence that the v p results, as estimated by MLP and NLLS, are not statistically different (p-value > 0.3), which effectively shows that both methods can be used interchangeably for blood volume fraction assessment. In the case of K trans , the difference in parameter estimates becomes statistically significant (p-value < 10 −3 ), although the medians remain relatively close with respect to the overall data variability. Estimates of T g and ∆ vary more significantly (cf. Figure 12c,d). The evidence against the null hypothesis of there being equal means in the MLP-and NLLS-based calculations is even stronger than for the K trans parameter (p-value < 10 −6 ). Apparently, the variation in the vascular impulse response function manifests itself in the CA concentration time courses in a way which is less evident in comparison to the other two PK model parameters. As such, the neural network could not precisely infer the relationship between the concentration time curves and the expected T g and ∆ values. This effect is consistent with the observed R 2 scores calculated on the test set upon MLP training. The VIRF timing factors are remarkably lower than the blood volume fraction and transfer constant. narrower for the neural network. Moreover, our ANNs ensure the better balance of the estimated values around the mean. The NLLS method, on the other hand, tends to overestimate the measured quantity, especially in the lower range of values. On the other hand, in the case of ANN, GFR was estimated without explicit determination of AIF for a given subject. It proves that the designed network was capable of extracting the general characteristics of the underlying perfusion process.

Discussion and Conclusions
The main goal of this study was to evaluate the possibility of estimating the parameters of the two-compartment filtration model using an artificial neural network. As shown in the Results section, especially in Table 4, overall the task was accomplished successfully. The estimated single-and twokidney estimates of glomerular filtration rate fall into the physiologically feasible range, and are close to the reference values calculated with the Trust Region Reflective algorithm (a non-linear leastsquares curve-fitting method), as well as those measured by blood tests. Comparison with the latter method, which serves as a gold standard in clinics, appears particularly optimistic. The mean difference between the corresponding measurements is to = 2.35 mL/min/1.73 m 2 , with the agreement interval = ± 36.16 mL/min/1.73 m 2 . This interval is even narrower than the one obtained for the NLLS method ( = ± 63.97 mL/min/1.73 m 2 ), which signals the greater precision of the designed MLP structures. This effect was achieved thanks to the usage of dropout layers between the hidden ones, and the training of the network with data vectors purposely corrupted with noise.
Although the trained MLP network performs relatively well in predicting K trans , vp and Δ parameters, it actually fails in estimating accurate values for the VIRF decay time constant Tg. The corresponding R 2 determination coefficients calculated for the test sets do not exceed the value of 0.64, which is remarkably less than for the other three parameters (0.85-0.92). Tg, however, does not affect the calculation of GFR, which depends only on K trans . In contrast with the NLLS procedure, a

Discussion and Conclusions
The main goal of this study was to evaluate the possibility of estimating the parameters of the two-compartment filtration model using an artificial neural network. As shown in the Results section, especially in Table 4, overall the task was accomplished successfully. The estimated singleand two-kidney estimates of glomerular filtration rate fall into the physiologically feasible range, and are close to the reference values calculated with the Trust Region Reflective algorithm (a non-linear least-squares curve-fitting method), as well as those measured by blood tests. Comparison with the latter method, which serves as a gold standard in clinics, appears particularly optimistic. The mean difference between the corresponding measurements is to µ d = 2.35 mL/min/1.73 m 2 , with the agreement interval κ = µ d ± 36.16 mL/min/1.73 m 2 . This interval is even narrower than the one obtained for the NLLS method (κ = µ d ± 63.97 mL/min/1.73 m 2 ), which signals the greater precision of the designed MLP structures. This effect was achieved thanks to the usage of dropout layers between the hidden ones, and the training of the network with data vectors purposely corrupted with noise.
Although the trained MLP network performs relatively well in predicting K trans , v p and ∆ parameters, it actually fails in estimating accurate values for the VIRF decay time constant T g . The corresponding R 2 determination coefficients calculated for the test sets do not exceed the value of 0.64, which is remarkably less than for the other three parameters (0.85-0.92). T g , however, does not affect the calculation of GFR, which depends only on K trans . In contrast with the NLLS procedure, a neural network does not actually fit a modeled curve to the input data. A failure in predicting a given parameter value does not necessarily cause an error in another parameter's estimation.
It is instructive to compare the results obtained in this study to similar measurements presented in [16], partly conducted on the same subjects. The scores presented therein were obtained using a conventional curve-fitting approach, and were divided with respect to examination session. The mean differences and limits of agreement reported in their study were µ d = 1.5 mL/min/1.73 m 2 and κ = µ d ± 43.2 mL/min/1.73 m 2 (Session 1), and µ d = 6.1 mL/min/1.73 m 2 and κ = µ d ± 31.9 mL/min/1.73 m 2 (Session 2). Hence, our MLP-based method estimates true GFR values at a comparable level, demonstrating either slightly better (in the case of Session 1) or poorer (Session 2) precision. Apart from different parameter estimation methods, the observed discrepancies may be caused by other factors, including various post-processing algorithms, e.g., image registration and segmentation. Above all, a different pharmacokinetic model has been utilized, namely the Sourbron's two-compartment separable model [7].
Regarding the PK models themselves, the most common of them are often claimed to be too simple to properly represent patient-specific biophysiological processes. Then, even the most precise parameter estimation may not lead to accurate GFR calculation. There have been several attempts to develop more advanced models, including multi-compartmental [32] or patient-specific models [33], as well as those with a modified functional form of CA retention in the renal tissue compartments [34]. As such, the reliability of our proposed MLP-based method for DCE-MRI processing should be further examined with respect to various mathematical models of kidney perfusion.
Another constraint that narrows the scope of the conclusions that can be derived from this study is that it presents methodological achievements tested solely on healthy subjects. Although we believe that the proposed methods of ANN-based PK model parameter estimation and AIF determination will manifest its applicability to diseased kidneys, this paper focuses exclusively on algorithms per se and numerical image processing, rather than clinical applications. Such an approach is common in numerous works published in the biomedical engineering domain, which postulate novel processing techniques, including those devoted to renal perfusion measurement [7,10,16,35]. The reasoning behind this approach states that the technical quality of the designed solution must be verified under well-defined and controlled conditions for a statistically representative data sample. Such conditions are ensured by a cohort of healthy volunteers. Besides, in the current study, the inclusion of abnormal cases would not necessarily provide additional insight into the behavior of the automatic determination of AIF. In the case of diseased kidneys, it should run in an unchanged manner, since tissue lesion is independent of the blood flow dynamics in abdominal aorta. On the other hand, taking into account renal impairments would involve experimenting with various pharmacokinetic models of perfusion, and it would remarkably extend the scope of this study.
To conclude, the performed experiments show the potential of neural networks as an alternative computational framework to the standard curve-fitting procedure in the context of quantitative perfusion estimation. The overall agreement with the gold standard method obtained by ANN was comparable to the results obtained using the non-linear least-squares approach. Moreover, GFR can be estimated without explicit determination of AIF for a given patient, which constitutes the major advantage of the proposed approach. Instead, AIFs must be found only for a limited number of studies included in the training data set. The described algorithm for AIF ROI determination allows one to determine it in an automatic manner, reducing the effect of the inflow artifact. It simultaneously increases the reliability and precision of the 2CFM model's parameter estimation. Hence, no additional process is required in assessing patient-specific hemodynamic parameters by utilizing either clinical examinations (MR angiography [36,37] or ultrasound [38]) or non-clinical methods (e.g., photopletysmographic arterial pulse wave measurement [39]) that would otherwise be necessary to properly trigger the DCE-MRI sequence.  Assuming the gradient echo sequence, C tissue in a given time step t is derived from the observed signal S(t) using the transformation C tissue (t) = R 1 (t) − 1/T 10 r 1 (A1) where T 10 is the longitudinal relaxation time before arrival of the tracer agent bolus, r 1 is the compartment-specific longitudinal relaxivity constant and In the equations above, α and TR are the MRI sequence parameters of flip angle and repetition time, S 0 is the signal baseline, i.e., before arrival of the CA bolus, and R 10 is the relaxation constant equal to 1/T 10 . Note that in the case of the tissue signal, the conversion function is valid only under the assumption of equal longitudinal relaxivities in the EEV and IV compartments.
where θ i,j ,θ i,j denote the actual and estimated value of a j-th parameter for an i-th data training vector, whereas M and N are the total numbers of parameters and data vectors, correspondingly. Since the model parameters have different scales, it was necessary to scale their values to the same range so as to balance their contribution to the overall loss function.