Pharmacokinetic (PK) models are used to extract physiological information from dynamic contrast enhanced magnetic reso- nance imaging (DCE-MRI) sequences. Some of the most common models employed in clinical practice, such as the standard Tofts model (STM) or the extended Tofts model (ETM), do not account for passive delivery of contrast agent (CA) through diffusion. In this work, we introduce a diffusive term based on the concept of effective diffusivity into a finite element (FE) implementation of the ETM formulation, obtaining a new formulation for the diffusion-corrected ETM (D-ETM). A gradient- based optimization algorithm is developed to characterize the vascular properties of the tumour from the CA concentration curves obtained from imaging clinical data. To test the potential of our approach, several theoretical distributions of CA concentration are generated on a benchmark problem and a real tumour geometry. The vascular properties used to generate these distributions are estimated from an inverse analysis based on both the ETM and the D-ETM approaches. The outcome of these analyses shows the limitations of the ETM to retrieve accurate parameters in the presence of diffusion. The ETM returns smoothed distributions of vascular properties, reaching unphysical values in some of them, while the D-ETM accurately depicted the heterogeneity of KTrans, ve and vp distributions (mean absolute relative difference (ARD) of 16%, 15% and 9%, respectively, for the real geometry case) keeping all their values within their physiological ranges, outperforming the ETM. Keywords Finite element method · Inverse analysis · Pharmacokinetic modelling · Magnetic resonance imaging 1 Introduction Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) has long been used as a clinical tool to study the vasculature of different tissues, especially tumours. Recently, it has been employed to assess the efficacy of antiangiogenic treatments in tumours [1 –3]. This technique consists on the acquisition of a series of T1-weighted images before, during, and after the intravenous injection of a con- trast agent (CA). As the CA reaches the tissue of interest via arterial inflow, it decreases the native T1 relaxation time, producing an increase in the measured signal intensity. The following removal of CA through venous blood flow results in a decrease in signal intensity, returning to its baseline value. For each voxel, the outcome is a signal intensity ver - sus time curve [4 ]. Three main approaches can be found in the literature to analyse these curves [5 ]: subjective analyses, based on the observation of DCE-MRI images by an experienced observer; semi-quantitative analyses, which extract param- eters from the signal-time curve, such as the initial slope, time to peak, etc; and finally, full-quantitative analyses. The last approach implies the conversion of signal intensity to CA concentration, which can be elaborated considering either linear or non-linear equations [6 ]. Once converted to CA concentration, pharmacokinetic (PK) models are applied to extract physiological information from the images. There are a wide range of PK models [7 ], considering different formulations and hypotheses to analyze the DCE data. One * José Manuel García- Aznar jmgaraz@unizar.es Diego Sainz-DeMena dsainz@unizar.es Wenfeng Ye wenfeng.ye@ansys.com María Ángeles Pérez angeles@unizar.es 1 Department of Mechanical Engineering, Aragon Institute for Engineering Research (I3A), University of Zaragoza, Zaragoza, Spain 2 ANSYS France, Villeurbanne, France 3850 Engineering with Computers (2022) 38:3849–3865 1 3 of the most common approaches is the so-called standard Tofts model [8 ] (STM). It is a compartmental model, which considers two different well-mixed compartments: intravas- cular and extravascular-extracellular space (EES), and the exchange of CA between them. In some highly vascularized tissues (such as liver tissue and some type of tumours), the voxel volume occupied by vessels is not negligible. Tofts et al. [9 ] formulated the extended Tofts model (ETM), which includes an additional variable that represents the contribu- tion of intravascular CA to signal enhancement. These models assume that CA can only reach the region of interest (ROI) through blood perfusion, neglecting passive delivery (spatial diffusion from adjacent regions, also known as inter-voxel diffusion, and convection due to gradients in interstitial fluid pressure). The hypothesis of well-mixed compartments assumes that CA is distributed uniformly within each compartment. This approximation also under - estimates the possible intra-voxel CA diffusion. Although previous works have studied the influence of this intra-voxel diffusion on fitted parameters [10, 11], this work will focus on inter-voxel diffusion. Even though convective effects may be relevant in some tissues, they will not be considered to reduce the complexity of the model. As several authors have pointed out [12– 17], assuming no inter-voxel CA diffusion can lead to errors in the quantities estimated by the model, especially in weakly vascularized zones, such as necrotic regions in a tumour. Previous works have developed different methods to include the diffusive phenomenon in the model. The models developed by Jia et al. [13] and Koh et al. [14] accounted for inter-voxel CA diffusion, but lacked information about other physiological properties, like those parameters defined in the STM or the ETM [8 , 9]. Pellerin et al. [15], first, and Fluckinger et al. [16], next, have created models that, keep- ing the formulation proposed by Tofts et al. [ 8], added a term to account for CA diffusion. Pellerin et al. [15] proposed a combined diffusion-perfusion (DP) that incorporated CA diffusion to the STM. However, since it required fitting all model parameters at once, the resulting computational cost was excessively high. Fluckinger et al. [16] solved this issue by implementing a voxel-wise approach to compute CA dif- fusion. To do so, they considered some simplifications, such as homogeneity in cellularity and diffusivity among neigh- bouring voxels, that might not apply to many types of tis- sues. Later, Cantrell et al. [17] proposed a diffusion compen- sated Tofts model and applied it to intracranial aneurysms. This approach considered a known diffusivity coefficient and handled separately the contributions due to diffusion and to extravasation to reduce the computational cost. Their method is based on the results obtained from the ETM, which may differ from true values in some tumours where CA diffusion plays an important role. In those cases, the applicability of the model is limited. More recently, Sinno et al. [18] have developed a new transport model (Cross Voxel Exchange Model, CVXM) that includes both the diffusion and the con- vection phenomena into the STM formulation. In their work, they quantified the error in the parameters retrieved by the STM due to ignoring passive transport mechanisms, show - ing the importance of convective and diffusive processes in DCE-MRI data. The implementation of the model was based on a discrete formulation, similar to previous works [15, 16]. They tested the model both on in silico and experimen- tal xenograft data, considering only 1D geometries. Except Cantrell et al. [17], other authors based their works on the STM, which is only accurate in poorly vascularized tissues [19]. To extend the range of application of the diffusion-cor - rected models, the present work is based on the ETM, which is accurate on both weakly and highly vascularized regions. Therefore, this study proposes a finite element (FE)-based model to include a diffusive term in the ETM [ 9] formula - tion with the aim of using this model to solve an inverse problem that provides an estimation of vascular properties of the tumour from DCE-MRI sequences. To achieve this aim, we first define the formulation of a diffusion-corrected extended Tofts model that accounts for the effect of CA dif - fusion and perfusion and we implement it on a FE model. Second, we develop a gradient-based optimization method to fit this model to the clinical imaging data to extract vascular properties. And third, we test the performance and accu- racy of this inverse method in different simulated theoretical cases. Hence, this approach aims to benefit from the com- putational efficiency of the finite element method (FEM) to fit all analysed volumes simultaneously while keeping an affordable computational cost. Thus, the modelling approach here proposed can be applied to analyze the heterogeneous behavior characteristic of tumours, avoiding simplifications and approximations that reduce the range of application of the model. 2 Materials and methods This section is divided into three different parts. In the first one, a whole description of the model formulation is pre- sented in detail. The second comprises the computational implementation of this model in the FE-based commercial software ANSYS (Ansys Inc., TX, USA). And finally, there is a complete description of the algorithm that solves the inverse problem of curve-fitting. 2.1 Diffusion‑corrected extended Tofts Model (D‑ETM) The transport mechanisms in biological tissues have long been of interest among physicians and researchers. Nicholson and collaborators [20– 25] studied thoroughly 3851 Engineering with Computers (2022) 38:3849–3865 1 3 this process, establishing a general formulation for the diffusion of particles in the brain [22]. In this work, they compared the densely packed cells of the brain and their interstitial spaces to a porous medium with two phases, an intra- and extracellular phase. They then formulated the equations of the diffusive process in brain based on the general equation of diffusion in porous media, which is a process that has also been widely studied by many authors [26– 28]. In the present work, this diffusion formulation is added to the general form of the ETM. There is a distinction between two different scales: a macro-scale, where the extravasation contribution to the concentration averaged in the representa- tive volume element (RVE) as well as the diffusion of CA between adjacent RVEs are studied (Fig. 1b); and the micro- scale, which is defined within the RVE and consists of a het- erogeneous distribution of cells, vessels and EES (Fig. 1a). The size of the RVE must be large enough to achieve length scale separation, containing sufficient number of cells and, at the same time, it must be small enough so that CA concen- tration can be averaged in it without adding significant error. Therefore, Eq. 1 describes the main transport mechanisms for a RVE: where: • Ct(x,t) is the total CA concentration, averaged in the RVE (VT ) From now on, we will refer to Ct(x,t) as Ct(x, t) to facili- tate the formulation • x is the coordinates vector in the macroscopic domain of the tissue • X is the coordinates vector in the microscopic domain of the RVE The relation between the macroscopic (x ) and the microscopic ( X) coordinate systems is defined through the homogenization process defined in Eqs. 2 and 3 • Deﬀ(x) is the effective diffusion coefficient for that RVE. Considering that each RVE contains a heterogeneous distribution of cells and vessels, we can compare the dif - fusive process of CA to the diffusion in a porous medium, as stated previously [22, 25–29]. A further detailed study of this variable is included below (Sect. 2.1.1) • KTrans(x) is the extravasation rate for that RVE(1)/u1D715Ct(x,t) /u1D715t=/u1D6C1⋅/parenleft.s2 Deﬀ(x)∇Ct(x,t)/parenright.s2 +KTrans(x)/parenleft.s2 Cp(t)−Ce(x,t)/parenright.s2 +vp(x)dCp(t) dt (2) Ct(x,t)=1 VT/uni222B.dspVTCt(X,x,t)dV • vp(x) is the volume fraction of blood plasma in the RVE • Cp(t) is the CA concentration in the blood plasma vol- ume, which follows the same temporal function in every RVE and only varies with time • Ce(x,t) is the averaged CA concentration in the extracel- lular subvolume (VE ) of the RVE. Considering that the microscopic scale is unknown (we have no information about the subvolume VE of the RVE), we assume the hypothesis of well-mixed compartments presented by Tofts [9 ]. This hypothesis considers an infinite diffusiv - ity coefficient of CA in the subvolume VE , what implies that there cannot be any spatial gradient of Ce in the sub- volume VE . Therefore: Fig. 1 a Insight on the different scales considered in the model. The RVE consists of a certain volume occupied by cells, another region which corresponds to the microvasculature and the rest of the volume (coloured in light yellow), which is the EES (VE ). The concentration of CA within the EES volume of the RVE is averaged (Eq. 3). The two different coordinate systems (x and X) are detailed. Panel b Dif- ferent processes taking place in the RVE (VT ) that are considered in the proposed model. Apart from the perfusion process (solid green arrows) and the contribution of microvasculature to the total con- centration, diffusion between adjacent RVEs (dotted blue arrows) is included 3852 Engineering with Computers (2022) 38:3849–3865 1 3 From now on, we will refer to Ce(x,t) as Ce(x,t) to facili - tate the understanding of the equations This well-mixed compartments hypothesis validates the Tofts equation for the compartmental model: where ve(x) is the volume fraction of the EES in the RVE, defined as ve=VE VT. Substituting Eq. 4 in Eq. 1, we obtain Eq. 5, which is the general form of the diffusion-corrected ETM (D-ETM), for - mulated in terms of the total concentration in the RVE: A table summarizing all model parameters, their definition and units can be found in Appendix A . 2.1.1 Effective diffusivity Different authors have formulated equations that relate the effective diffusivity to different geometrical characteristics of the porous material [29– 31]. Given that the geometrical structure of the solid phase (cells and vessels) is unknown, we assume that the equivalent diffusivity is related to the tortuos - ity [22]. Tortuosity ( /u1D706 ) quantifies the increase in path-length of a diffusing particle due to the existence of obstacles in its way [ 32]. It is defined as: Tao and Nicholson [ 25] used the Monte Carlo method to simulate the diffusion of particles on different porous struc - tures and different cell geometries and fitted the simulated data to obtain the value of Deﬀ for each case. Applying Eq. 6 to these values, they obtained a value of tortuosity ( /u1D706 ) for each simulated structure and cell geometry and found that tortuosity was independent on the considered cell geometry and was only dependent on the porosity ( /u1D700 ) of the structure: Substituting Eq. 7 in Eq. 6 we obtain the definition of an equivalent diffusivity (Deﬀ ). Considering that porosity ( /u1D700 ) is defined as the volume fraction of ”empty” space in the material, its equivalence to ve is straightforward:(3) Ce(x,t)=1 VE/uni222B.dspVECe(X,x,t)dV (4) Ct(x,t)=ve(x)Ce(x,t)+vp(x)Cp(t) (5)/u1D715Ct(x,t) /u1D715t=/u1D6C1⋅/parenleft.s1Deﬀ(x)∇Ct(x,t)/parenright.s1 +KTrans ve(x)/parenleft.s1Cp(t)(ve(x)+vp(x)) − Ct(x,t)/parenright.s1+vp(x)dCp(t) dt (6) /u1D706=/uni221A.s2 D Deﬀ (7) /u1D706=/uni221A.s2 3−E 2Where D is the diffusion coefficient of CA in free medium, which is known. Using this formulation, we can consider different effective diffusion coefficients per element without adding more unknowns to the model. 2.2 Forward FE model Once the theoretical formulation is defined, the next step is to implement it into the FE software to generate CA vs time curves from given sets of parameters. Equation 5 is imple- mented into ANSYS diffusion module [33], including the extravasation term and the contribution of the blood plasma fraction as non-linear generation terms. The numerical formulation obtained for Eq. 5 is defined as: Where, for an element e: • ce is the nodal concentration vector and ̇ce is its temporal derivative • n are the element shape functions • Cd e is the element diffusion damping matrix and it is defined as: • Kd e is the element diffusion conductivity matrix and it is defined as: where D is the diffusion coefficient matrix, further defined • rg e is the element diffusing substance generation load vec- tor and it is defined as: where g is the generation load vector. The integration volume Ωe corresponds to the volume of the finite ele - ment e. The time discretization scheme is derived from the backward Euler implicit method and is defined as: Where:(8) Deﬀ=2D 3−ve (9) Cd ėce+Kd ece=rg e (10) Cd e=/uni222B.dspΩennTd(Ωe) (11) Kd e=/uni222B.dspΩe/parenleft.s1∇nT/parenright.s1TD∇nTd(Ωe) (12) rg e=/uni222B.dspΩegnTd(Ωe) (13) un+1=un+/u1D703Δt ̇un+1+(1−/u1D703)Δt ̇un 3853 Engineering with Computers (2022) 38:3849–3865 1 3 • un is the nodal degree of freedom (DOF) values at time tn • /u1D703 is the transient integration parameter. If /u1D703=0 an explicit algorithm is used, whereas if /u1D703=1 an implicit algorithm is employed • Δt=tn+1−tn is the time step size • ̇un is the time rate of the nodal DOF values at time tn , computed at previous time point Implementing this time discretization scheme on Eq. 9 leads to Eq. 14: Where R is the substance generation matrix and Kd is the assembled diffusion conductivity matrix. The terms ct+1 and ct are the nodal concentration vectors resulting from evaluat- ing the nodal concentration matrix ( C ) at given time points. kTrans , ve and vp are the nodal variables vectors. Finally, ct+1 p and ̇cpt+1 are the scalars resulting of evaluating the arterial input function (AIF) and the derivative of the AIF vectors ( cp and ̇cp ) at given time points. To distinguish between the regular dot product and the Hadamard operations, we employ a specific notation for the pointwise product ( ⊙ ) and division ( ⊘ ). The diffusion coefficients matrix D for each element is: Where veelem is the element average of the nodal ve vari - able. Considering as a first approach isotropic diffusion, Dxx=Dyy=Dzz=D. Due to the dynamic nature of the physical process, tran- sient effects were included. Although the transient integra- tion parameter was set to /u1D703=1 , which makes the solution unconditionally stable, the influence of the transient effects still affected the accuracy of the solution depending on the time step ( Δt ). We conducted several experimental error analyses and concluded that a time step of 1s ensured the consistency of the solution, while keeping an affordable computational cost. Initial values were provided to initialize the transient sim- ulations. Both the nodal concentration and its time derivative were set to zero for t=0s ( C(t0)= ̇C(t0)=0). The AIF is interpolated to match this time resolution. Since both generated curves and fitted curves are obtained using the same model, no interpolation is needed to perform the fitting process.(14)Cdct+1−ct Δt=−Kdct+1+R/parenleft.s2 (kTrans⊘ve)⊙(ve+vp)ct+1 p/parenright.s2 −R/parenleft.s1(kTrans⊘ve)⊙ct+1/parenright.s1+Rvṗcpt+1 ⎡ ⎢ ⎢ ⎢ ⎢⎣2Dxx 3−velem e00 02Dyy 3−velem e0 002Dzz 3−velem e⎤ ⎥ ⎥ ⎥ ⎥⎦Regarding the boundary conditions, we assumed that no CA could diffuse across the boundaries of the geometry, considering the tumour as an ”isolated” entity with respect to adjacent tissues. Although it may not be biologically cor - rect for tissues that are not surrounded by physical barriers, this condition facilitates the formulation and is consistent with literature [8 , 9, 15, 16]. 2.3 Solving the inverse problem Imaging data is usually processed as a voxelized geometry. However, the method proposed here can be applied to any complex geometry, like those obtained from imaging seg- mentation. Although the conversion from a voxelized discre- tization to a FE mesh approximates the data, since the spatial resolution of the latter is expected to be equal to or higher than the spatial resolution of the former, the error added to the data is small. After implementing the D-ETM equation in ANSYS, we can simulate the CA transport for a certain tissue. There- fore, the next step is fitting the model to the concentration- time curves obtained from imaging. This inverse problem of curve fitting is solved using an iterative method based on the non-linear least squares method [34]. Being a gradient based method, it needs to compute the derivative of the total concentration C (t) with respect to each nodal parameter, known as the Jacobian matrix. Typically, this matrix is obtained numerically using the finite differences method. This option leads to an excessively high computational cost, since the number of simulations ( Totalsim ) needed to compute this matrix increases linearly with the number of time points and nodal variables and quadratically with the number of nodes. where Nnodes is the number of nodes in the model, Ntime is the number of time points and Nparameters is the number of unknowns to be fitted per node. Even for small cases, this numerical analysis requires an excessive computational cost that cannot be considered. Therefore, we propose an alternative semi-analytical computation of the Jacobian matrix. We can re-write Eq. 14 as: where:(15) Totalsim=Nnodes2NtimeNparameters (16)Act+1=Δ tR/parenleft.s2 (kTrans⊘ve)⊙(ve+vp)ct+1 p+vṗcpt+1/parenright.s2 +Cdct (17) A=Cd+Δ tKd+Δ tR(kTrans⊘ve) 3854 Engineering with Computers (2022) 38:3849–3865 1 3 The derivative of C with respect to the different parameters can be obtained from Eq. 16 by applying the product rule on the left hand side of the equation: • For KTrans : • For ve : where Kd(ve)=/u1D715Kd /u1D715ve . Kd is the assembled matrix of ele- mental diffusion conductivity matrices, which are defined by Eq. 11. To build these matrices, we define D as: Where, for an element e, vee=1 N∑N i=0vei . N is the num - ber of nodes in the element and vei are the nodal values of ve . The derivative of Kd e with respect to nodal ve ( Kd e(ve) ) is exported from ANSYS by deriving the components of D with respect to the nodal ve values, obtaining: Each elemental matrix Kd e(ve) is generated introducing this matrix in Eq. 11. Finally, these elemental matrices are assembled to build the global matrix Kd(ve). • For vp : To compute Eqs. 18, 19, 20, matrices Cd and R are exported from ANSYS along with Kd(ve) . After solving these equa- tions, the Jacobian matrix is obtained by concatenating the matrices computed in Eqs. 18, 19, 20:(18)A𝜕ct+1 𝜕kTrans=Δ tR⊙/parenleft.s2 ct+1 p(ve+vp)⊘ve−ct+1⊘ve/parenright.s2 +Cd𝜕ct 𝜕kTrans (19)A𝜕ct+1 𝜕ve= −Δt Kd(ve)⊙ct+1 +Δ t/parenleft.s2 R⊙(kTrans⊘ve2)⊙ct+1−vpct+1 p/parenright.s2 +Cd𝜕ct 𝜕ve ⎡ ⎢ ⎢ ⎢⎣2D (3−vee)00 02D (3−vee)0 002D (3−vee)⎤ ⎥ ⎥ ⎥⎦ ⎡ ⎢ ⎢ ⎢ ⎢⎣−2D N(3−1 N∑N i=0vei)200 0−2D N(3−1 N∑N i=0vei)20 00−2D N(3−1 N∑N i=0vei)2⎤ ⎥ ⎥ ⎥ ⎥⎦ (20)A𝜕ct+1 𝜕vp=Δ tR⊙/parenleft.s2 (kTrans⊘ve)ct+1 p/parenright.s2 +Rİcpt+1+Cd𝜕ct 𝜕vpThis matrix is included in the least-squares method to fit the CA concentration-time curves to Eq. 5. We tested this analytical computation in small-sized models, comparing it to the numerically computed Jacobian. Results were promis- ing, as both matrix and convergence were similar between the two methods. Further information about these results can be found in Appendix B. The optimization algorithm developed takes as inputs an initial set of values of the unknowns, as well as the curves to be fitted. The model proposed defines a set of three param- eters (KTrans , ve and vp ) per node. This implies that the algo- rithm needs to fit three times the number of nodes simulta- neously. Starting from the initial seed provided, a forward simulation of the D-ETM is executed, obtaining a first set of curves that are then used to initialize the cost function. Then, the numerical matrices needed for the Jacobian computation ( Cd , R and Kd(ve) ) are exported by running several scripts on ANSYS. Once the Jacobian is obtained, the minimization solver computes the updated set of parameters, finishing the first iteration. A schematic pseudocode of this process is presented in algorithm 1. The cost function (CF) defined for the method is a stand- ard sum of squared differences (Eq. 21), which has proven to be effective in the simulated cases: Where Ntime is the total number of time points, yi is the value of the reference concentration at time point i and fi(xj) is the value at time point i of the curve obtained from running a simulation of D-ETM with given xj vector of parameters. Other cost functions that only consider a fraction of the time points, like those proposed in [15], were tested. Never - theless, they showed convergence problems, leading to the algorithm getting caught in local minima distant from the global optimum. The solver method chosen to perform the minimization was the Trust Region Reflective (TRF) algorithm [35]. Although the commonly used Levenberg-Marquardt (LM) [36] is also suitable for our purposes, the TRF algorithm (in its Scipy [37] implementation) handles sparse matrix. Given the nature of our problem, working with matrices in sparse format for the Jacobian computation dramatically reduced the use of system memory (around 98% reduction). We compared the performance of the LM algorithm against J=⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣/u1D715c1 1 /u1D715kTrans 1…/u1D715c1 nnode /u1D715kTrans nnode/u1D715c1 1 /u1D715ve1…/u1D715c1 nnode /u1D715vennode/u1D715c1 1 /u1D715vp1…/u1D715c1 nnode /u1D715vpnnode ⋮ /u1D715cntime 1 /u1D715kTrans 1…/u1D715cntime nnode /u1D715kTrans nnode/u1D715cntime 1 /u1D715ve1…/u1D715cntime nnode /u1D715vennode/u1D715cntime 1 /u1D715vp1…/u1D715cntime nnode /u1D715vpnnode⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ (21) CF=Ntime/uni2211.s1 i=0(yi−fi(xj))2 3855 Engineering with Computers (2022) 38:3849–3865 1 3 the TRF (converting the matrices to dense format for the LM method). Both methods retrieved accurately the reference values and the number of iterations needed was similar in both cases. The proposed curve generation and curve fitting pro- cesses were executed in a cluster composed of 480 CPUs and 1088 GB of RAM. Specifically, the resources employed for the benchmark problem and the real tumour geometry were: 16 CPUs and 16GB RAM and 24 CPUs and 32GB RAM, respectively. The optimization time of the benchmark problem was 1.5h on average (fitting 360 time points on 339 nodes) and for the tumour geometry was around 10h on average (fitting 360 time points on 955 nodes). The optimi- zation algorithm was developed in Python and APDL and the forward simulations were generated on ANSYS 2019R2. orders of magnitude) to parallelize several 2D simulations rather than studying all of them as a 3D case. The geometries correspondent to both cases were meshed using linear quadrilateral diffusion elements (PLANE238 in ANSYS manuals ) of size 0.15 mm, which generates a mesh of higher density than of most of the clinical and experimen- tal imaging data. Similarly to previous works [10], the diffusion coefficient D was given a value of 2.6E-04 mm2/s, which is consistent with experimental measurements [14, 40]. Considering the formulation of the diffusive term of the model proposed here (Eq. 5), the diffusion coefficient was constant in both cases. The AIF used in the simulations (Fig. 2) was inspired on the clinical data collected within the European research project PRIMAGE [38]. Fig. 2 AIF used for the simulations 35 Algorith m1 Pseudo codeoftheoptimizatio nprocess 1:Generate aninitia lsetofparameters (P0) Initia lseed 2:while (costfunctio n>Threshold 1)or(δ>Threshold 2)do Convergenc ecriteria. 3: RunforwardFEmodelwiththecurren tsetofparameters (Pi) Obtain ﬁttedcurves 4: Evaluate CF 5: Exportnumerica lmatrices from ANSY SExportCd,RandKd(ve) 6: Comput eJacobia nmatrix 7: Deﬁne newsetofparameters Pi+1=Pi+ where δiscomputed by theTRF[]viathegradie nt. 8:endwhile 3 In silico simulations To test the performance of the optimization algorithm pro- posed and to compare the D-ETM to the ETM, two different sets of simulated tissue concentration time courses have been generated from the forward FE model of the D-ETM and some sets of known parameters (KTrans , ve , vp and D). The first case corresponds to a simple two-dimensional (2D) geometry, where the distribution of parameters gener - ates a CA distribution that is completely dependent on diffu- sion. The second case, on the other hand, is based on a more complex geometry inspired in clinical data of tumours [38]. Although the proposed method can be applied to both 2D and three-dimensional (3D) cases, given the common resolution of imaging data (slice thickness is usually several times the pixel size), we focused on 2D geometries to keep the computational cost affordable. Due to this size differ - ence, the effects of out-of-plane diffusion are expected to be negligible in comparison with the in-plane diffusion [16]. In terms of computation, it is dramatically faster (up to two 3856 Engineering with Computers (2022) 38:3849–3865 1 3 3.1 Benchmark problem This case is based on the one proposed by Pellerin et al. [15], which was also included in the work of Fluckiger et al. [16]. Given the importance of these contributions into the study of CA diffusion process, we consider this case as a benchmark problem for PK models that incorporate diffusion. Hence, we simulated a slice of a circular, radially sym- metric, tissue. The definition of both the geometry and the distribution of parameters aims to generate an extreme case of a distribution of CA that is diffusion-limited. To do so, the circle is divided into two different regions: a highly perfused rim and a necrotic core. Initially, the parameters values chosen were similar to those used by Fluckiger et al. [ 16]: KTrans=0.2 min−1 in the rim and 0.05 min−1 in the core; a constant value of ve equal to 0.5 in the whole model; and, finally vp=0.05 in the rim and 0.005 in the core. When running the optimization algorithm on this case, we observed that, although the cost function was reduced to values close to zero, the parameters returned were dif- ferent from the true ones. Upon closer examination, we noticed the existence of several local minima close to the global optimum. Figure 3 shows two different nodes where the true values of the model parameters were the same, and so were the generated CA curves. However, although the fitted curve is almost identical to the reference one, the two sets of fitted parameters are different between them and both differ from the true values. It must be noted that some of the parameters retrieved are unphysical: KTrans and vp below zero in the second case of Fig. 3. Even when applying bounds to keep the parameters within physical ranges ( KTrans greater than zero and ve and vp between zero and one), the optimization algorithm still got caught in a local minimum. This meant that the success of the optimization process was dependent on the initial seed. The curves shown in Fig. 3 were obtained using as initial seed a set of values that was a random distribution of values between 0.4 min−1 and 0, for KTrans ; between 0.2 and 0.8 for ve ; and between 0 and 0.1 for vp . In this case, the optimization method retrieved accurate results only if the initial seed was very close to the true values. We attribute this problem to the numerical instability pro- duced by the radial symmetry, both in geometry and param- eters, which produces a set of identical curves at nodes with similar values. To prove this hypothesis, we created two additional simulations. In the first one, the parameters were kept constant for each region, while the geometry was modi- fied to suppress the axial symmetry. The results obtained from these simulations are detailed in Appendix C. In the second simulation, the geometry was preserved, while the parameters were changed by random distributions of values within a range. Such, KTrans was assigned values between 0.25 min−1 and 0.15 min−1 in the rim and between 0.05 min−1 and 0 min−1 in the core. vp , on the other hand, was given val- ues between 0.08 and 0.03 in the rim, and between 0.01 and 0 in the core. ve maintained its original values. The initial seed was the same as in the previous simulation. Results obtained with the ETM on this simulation are consistent with previous works [15, 16]. The KTrans , mainly, and the vp maps, to a lesser extent, show an averaging pat- tern with respect to the reference values. Besides, the ve map returns unphysical values (ve>1 and ve<0) in the necrotic region. Even if the KTrans and vp in the necrotic region are not exactly zero, the enhancement curve of these nodes is completely dependent on diffusion. Thus, these curves can- not be accurately fitted to the Tofts formulation. Besides, the Fig. 3 Two of the sets of curves (simulated and fitted) correspondent to some of the nodes within the necrotic region of the first simula- tion of the benchmark problem. The true values for both nodes are: KTrans=0.05 min−1 , ve=0.5 and vp=0.005. The fitted values are: in the first case, KTrans=0.15 min−1 , ve=0.41 and vp=0.01; and, in the second case: KTrans=-0.09 min−1 , ve=0.25 and vp=-0.001. These results dem- onstrate the convergence of the inverse method to a local minima 3857 Engineering with Computers (2022) 38:3849–3865 1 3 vascularized regions adjacent to the necrotic ones are also influenced by these values, since there is diffusion of CA from the former to the latter. This diffusive process results in a lower CA concentration, which is then fitted to KTrans values below the true ones. vp is not as affected by this effect as KTrans is, resulting in a better fit for this parameter. The distributions obtained with the D-ETM are very accurate for KTrans and vp , while the fitted ve map shows higher error, especially in the necrotic region. This error is related to the influence of the variable in the global solution. If we take a look at Eq. 19, we can see that the value of this derivative is dependent on the value of KTrans and vp . Therefore, in those necrotic regions, where these parameters present values close to zero, the derivative depends only on the diffusive term. Due to the definition of this term (Eq. 6), the influence of ve is limited. Thus, because of its low effect on the global solution, the optimization algorithm is not able to retrieve accurate results for ve , especially on necrotic regions. A quantitative comparison between the outcome of both models is presented in Table 1. Since both the true and initial values for the simulations are a function of randomness, ten cases with different true and initial sets of parameters were tested to validate the robustness of the method. Considering that reference values for KTrans and vp reached zero in the necrotic region, the use of relative error metrics is unfeasi- ble. The absolute error, measured as the absolute difference between the fitted and the reference value, was selected to compare the performance of both models. Different error thresholds were defined to compare the performance of both models. Threshold for KTrans was set to 0.01 min−1 , which is the maximum precision of the DP model [15]. Similarly, the threshold for vp was set at 0.001, a value that can be considered as sufficient precision for this kind of models. Due to the impact of ve maps on the global solution, its threshold was set to 0.15. The D-ETM clearly outperforms the ETM, especially on KTrans maps. While only 19% of nodes fitted by the ETM are within a 0.01 min-1 range from the reference value, 72% of those retrieved by the D-ETM fall into that range. Besides, the KTrans mean absolute error in the ETM is around four times higher than the one correspondent to the D-ETM.Due to the great effect of the unphysiological values of ve returned by the ETM on the absolute error, this metric will not be used to compare the models. The fraction of nodes whose error is below 0.15, however, is not affected by these values. While only 55% of the values obtained using the ETM are within the error range, almost 78% of the ones retrieved by the D-ETM fall into this range. Despite experiencing difficulties retrieving the correct ve maps, the D-ETM shows a great improvement with respect to the ETM. Besides, all of the ve values fitted using the D-ETM were within the physiological range [0,1]. The mean absolute error of D-ETM vp maps is around half the error obtained by the ETM. Nevertheless, this parameter does not seem to be as affected by diffusion as the other two. Although the error retrieved by suppressing the homoge - neity in the distributions of parameters (Fig. 4) was lower than the error obtained by removing the axial symmetry in geometry (Appendix C), these results demonstrate that the combination of both factors was causing the convergence of the algorithm to local minima. 3.1.1 Analysis of the mesh effect The influence of mesh size on the convergence of both the forward and the inverse models was tested using the bench- mark problem geometry. Two different meshes were gener - ated: the first one discretized the geometry using 0.15 mm size elements, while the element size on the second one was half that value. The number of nodes on the two simulations were 339 and 777, respectively. The same type of elements (linear quadrilateral diffusion elements) was used on both models. On both cases, the sets of true and initial values were random distributions between the ranges defined previ- ously. Just as in the benchmark problem, ten different simu- lations were generated for the finer mesh, to eliminate the influence of randomness on the result. The results presented in table 2 show that the method reduces the error when refining the mesh. Nevertheless, the slight increase in accuracy does not justify the greater computational cost associated to the finer mesh. The finer Table 1 Comparison of error metrics between the D-ETM and the ETM for the benchmark problem KTrans absolute error is measured in min−1 . Results were computed from ten simulationsD-ETM ETM Absolute error Fraction of nodes error < thresholdAbsolute error Fraction of nodes error < threshold Mean SD Mean (%) SD (%) Mean SD Mean (%) SD (%) KTrans 8.14E-03 1.18E-03 71.54 6.03 3.71E-02 6.93E-04 19.03 2.19 ve 9.53E-02 4.42E-03 77.58 1.42 7.78E+06 5.46E+05 54.57 0.53 vp 3.79E-04 7.45E-05 94.13 3.02 7.08E-03 2.71E-04 74.93 1.35 3858 Engineering with Computers (2022) 38:3849–3865 1 3 mesh model needed twice the time of the original model to fit the curves. Therefore, the element size selected for the simulations was 0.15 mm, a tradeoff between accuracy and computational cost. 3.2 Real tumour geometry A second set of simulations was generated to test the per - formance of the D-ETM in real geometries and vascular properties with a heterogeneous distribution. The geom- etry corresponded to a tumour slice of around 20 mm2 , while the vascular properties were inspired by clinical data [38]. The vascular properties distribution is divided into three different zones: a highly perfused region, an intermediate region and a necrotic region. Depending on the zone, the assigned parameters were: for KTrans , random values between 0.4 and 0.3 min−1 , between 0.25 and 0.1 min−1 and between 0.05 and 0 min−1 , for the three respective regions. ve values were randomly selected within a range between 0.85 and 0.75 for the necrotic region and 0.6 and 0.4 for the other two. Finally, vp random distribution ranged from 0.08 to 0.03 for the highly perfused nodes, between 0.05 and 0.03 in the intermediate region and between 0.01 and 0.005 in the necrotic one. The convergence of the inverse method was tested by pro- viding random distributions of parameters as initial values. The ranges for KTrans , ve and vp were the same as in the previ - ous case: [0,0.4] min−1 , [0.2,0.8] and [0,0.1], respectively. Following the procedure described in the previous case, ten simulations with different reference and initial values were generated. The heterogeneous reference maps generated clearly expose the limitations of the ETM in accurately captur - ing the KTrans distributions (Fig. 5). While the D-ETM pro- vides an almost exact distribution for KTrans and vp and an Fig. 4 Second simulation of the benchmark problem. Com- parison between the reference values and the parameters returned by the D-ETM and the ETM. Results show that the D-ETM accurate retrieves the distribution of KTrans and vp , while the ETM show an averaging pattern, especially for KTrans . Although not as accurate as the other parameters, the v e map returned by the D-ETM is within the physiological range [0,1], while the distribution obtained from the ETM reaches values close to infinity in the necrotic coreKTrans ve vpTrue V alues D-ETM ETM Table 2 Comparison of error metrics between two different meshes KTrans absolute error is measured in min−1 . Results were computed from ten simulations0.15 mm element 0.075 mm element Absolute error Fraction of nodes error < thresholdAbsolute error Fraction of nodes error < threshold Mean SD Mean (%) SD (%) Mean SD Mean (%) SD (%) KTrans 8.14E-03 1.18E-03 71.54 6.03 7.65E-03 7.43E-04 72.5 4.12 ve 9.53E-02 4.42E-03 77.58 1.42 9.11E-02 3.96E-03 80.35 1.41 vp 3.79E-04 7.45E-05 94.13 3.02 3.24E-04 3.78E-05 97.03 1.85 3859 Engineering with Computers (2022) 38:3849–3865 1 3 acceptable ve map, the ETM tends to homogenize the KTrans , failing to depict the highly perfused regions, as well as the necrotic ones. This effect is particularly visible in those zones where two of these regions are adjacent. The error metric employed in this case was the abso- lute relative difference (ARD), calculated between true and fitted parameters. Table 3 shows the metrics correspond- ent to the average of ten different simulations of the real geometry case. In this case, where necrotic zones were not as large as in the previous case, the ETM accurately retrieves the ve map (mean ARD is 11% and 86% percent of nodes have an ARD below 20%), depicting the increase in ve correspondent to these necrotic regions (Fig. 5). The D-ETM, on the other hand, gives a good fit for the ve map (mean ARD is 16% and the ARD is below 20% in 77% of the nodes), although it is not as accurate as the ETM in fitting the ve values in necrotic nodes (Fig. 5). KTrans and vp distributions obtained from the D-ETM are almost an exact representation of the reference maps (Fig. 5), as evi- denced by the metrics obtained (Table 3). The mean ARD is 16% and 9% for KTrans and vp , respectively. The ETM, nonetheless, show higher error for these two parameters (mean ARD of 148% and 195% for KTrans and vp , corre - spondingly, and for both maps only 40% of nodes have an ARD below 20%). In their work, Pellerin et al. [15] tested the performance of their model in a simulated case similar to the real geom- etry case presented in here. The DP model obtained a mean ARD of 16% for KTrans and 17% for ve , with 73% of the KTrans values and 77% of the ve values within 20% of the true values. The D-ETM has obtained identical values for ve and a similar mean ARD for KTrans , improving the fraction of values whose ARD is below 20%. 3.2.1 Influence of noise To test the robustness of the D-ETM to the addition of noise, several simulations were conducted. Starting from a set of reference values similar to those generated in the last case, experimental levels of noise were added to the Fig. 5 Real tumour geometry. Reference values and results of the D-ETM and the ETM for each of the parameters. Result show that the D-ETM accurate captures the heterogeneity of the distribution of parameters, while the ETM tends to aver - age the values. The maps of vp are the least sensitive to this phenomenon. The v e distribu- tion obtained from the D-ETM shows a more accurate fit than the one obtained in the bench- mark case KTrans ve vpTrue V alues D-ETM ETM Table 3 Error metrics comparison between the D-ETM and the ETM for the second case Results were computed from ten simulationsD-ETM ETM ARD Fraction of nodes ARD < 20%ARD Fraction of nodes ARD < 20% Mean (%) SD (%) Mean (%) SD (%) Mean (%) SD (%) Mean (%) SD (%) KTrans 16.37 3.67 87.21 2.67 148.42 47.37 39.85 1.6 ve 15.54 1.02 76.87 2.69 11.36 0.22 85.71 1.07 vp 8.59 4.45 95.56 3.02 194.85 101.26 40.13 1.53 3860 Engineering with Computers (2022) 38:3849–3865 1 3 generated curves. These levels were defined using a gauss- ian distribution with an standard deviation (SD) equal to a fraction (1%, 2.5% and 5%) of the highest concentration reached in the curve. The results of these simulations presented in Fig. 6 show the effect of noise on both the D-ETM and the ETM. Despite showing higher error for noise-free simulations, the parameters returned by the ETM seem to be unaffected by noise, since the histograms show almost no difference in the distribution of ARD between the cases with different levels of noise. The added noise shows greater influence on the parameters obtained fitting the D-ETM. Noise seems to have the greater effect on vp and the slighter effect on ve . KTrans , for its part, shows a slight disturbance when the D-ETME TM Fig. 6 Influence of noise on the accuracy of the D-ETM (left) and the ETM (right). Results show that noise has greater influence on the D-ETM, particularly on vp . Although the ETM seems to be unaf-fected by noise, for noise values below or equal to 2.5% the D-ETM still performs better. Even for the maximum levels of noise (5%), the error in the D-ETM solution is similar to the one in the ETM 3861 Engineering with Computers (2022) 38:3849–3865 1 3 noise level is low (1% and 2.5%), having more than 70% of its values within a 20% range from the true values. When the added noise reaches the maximum value, this percent- age drops dramatically to 30%. These results are consistent with those obtained by Pel- lerin et al. [ 15]. In their work, the authors attributed this different effect of noise on parameters depending on the influence of each parameter on the different parts of the CA curve. Thus, ve regulates the last part of the curve, where the concentration is close to the maximum and, therefore, it is less sensitive to noise. KTrans , and even more vp , affect the initial part of the curve. Consequently, the influence of added noise is greater on these parameters. The great number of variables to be fitted simultane- ously makes this model, as well as the model developed by Pellerin et al. [15], more sensitive to noise. Therefore, when fitting the model to experimental data it must be ensured that the signal-to-noise ratio (SNR) of the data meets the model requirements. 4 Discussion The use of DCE-MRI sequences to assess the efficacy of antiangiogenic therapies in tumours [1 –3] increases the need for PK models that retrieve vascular properties as accurate as possible. Several authors have pointed out the limitations of the widely used standard and extended Tofts models when the CA reaches the region of interest (ROI) through passive delivery [12– 18]. These models return an inaccurate estimation of KTrans as well as unphysical values for ve in those regions within the ROI where the active delivery of CA is low or non-existent (necrotic zones). Dif- ferent models and methods have been developed to assess the effect of diffusion and to develop PK that accounted for this process [15– 18]. These works exposed the men- tioned limitations and proposed different approaches to include the effects of diffusion into the STM and the ETM. Pellerin et al. [15] were the first to include a diffusive term into the STM. The DP model proposed showed an improvement in parameter accuracy with respect to the original STM in those regions where passive delivery of CA was significant. One of its major limitations was the high computational cost associated to the simulated annealing algorithm, since the model had to fit all voxels simultaneously. Fluckiger et al [16] added some hypoth- eses to the DP model (homogeneity in the diffusive coef- ficient and ve between neighbouring voxels) that allowed them to compute the effects of diffusion while fitting each voxel separately and, therefore, reducing substantially the computational cost. However, this homogeneity hypothesis may not be suitable for some kind of tumours (such as neuroblastoma, which is a type of cancer characterized by its high heterogeneity [41]). Cantrell et al. [17] based their formulation of the diffusive term on the one proposed by Pellerin et al. [15] and proposed a diffusion-compensated Tofts model (DC-Tofts). Their work consisted on obtain- ing the ETM parameters from CA curves, then computing the diffusion contribution from these fitted parameters to generate a new set of curves and, finally, fit again these last curves with the ETM. Although this model was suitable for intracranial aneurysms, the approach used may not be appropriate for other type of lesions where necrotic zones are present, since the initial fit of the ETM would retrieve unphysical values of ve that would condition the following diffusion computation. All these works follow the same formulation of the diffusive term (revised thoroughly in [15]). This formulation assumes that CA can diffuse freely through each of the voxel faces (this means that no obsta- cles, such as cells, are present in that faces). Depending on the cellularity level and the voxel size, this assumption may not be valid. The latest contribution by Sinno et al. [18] proposes a modified formulation of the diffusive term that, as well as the formulation presented in this work, avoids this simplification. In their work, the diffusivity coefficient is considered unknown but constant through the domain. As stated previously, this hypothesis may not be valid for some tissues, especially in tumour tissues, which are characterized by their heterogeneity. The formulation of the diffusion process here proposed avoids this assumption by embracing the concept of effec - tive diffusivity [22]. This concept implies that the diffusion of agents within biological tissues is similar to the diffusion of an agent in a porous medium [26– 28]. Apart from pro- viding a more accurate description of the diffusive process, this hypothesis links the effective diffusivity to the volume fraction of the extravascular-extracellular space (ve ), avoid- ing the generation of a new parameter (D) that needs to be fitted or extracted from data or literature, as it is the case in previous works [14– 18]. The performance of the D-ETM was tested using two different in silico cases. The first one, a benchmark problem derived from literature [15, 16], exposed the limitations of the ETM in those regions were passive delivery of CA is the main transport mechanism. Results showed the improved accuracy of the model parameters returned by the D-ETM, which were very close to the reference values (mean abso- lute error for K Trans , ve and vp were 0.008 min−1 , 0.095 and 0.0004, respectively). A second case inspired in real tumour geometries and parameters was analysed. Again, the ETM performed poorly, returning an incorrect estimation of KTrans . In both simulations, the KTrans distribution obtained by the ETM appeared averaged with respect to the reference maps, adding significant error to the parameters. D-ETM param - eters were almost identical to the true values, accurately 3862 Engineering with Computers (2022) 38:3849–3865 1 3 depicting the heterogeneous distribution of values. Besides, Fig. 6 shows that the error obtained by the D-ETM follows a Pareto-like distribution for noise levels below 2.5%, what means that the smallest errors are obtained for most of nodes, while the ETM does not follow this pattern (except for ve ), even for noise-free simulations. The results obtained from the D-ETM in both cases show lower accuracy in ve maps with respect to the other vari- ables, especially on necrotic regions, where KTrans values are very close to zero. This is due to the modest influence of ve on the global solution in these regions. To solve this issue, future works could develop an alternative expression for the effective diffusivity where ve would have more influence. Nevertheless, the accuracy of the ve maps obtained from the D-ETM is acceptable, keeping all of their values within the physiological range. Moreover, the ve maps retrieved by the D-ETM on the real tumour geometry show an accuracy similar to the DP model [15]. Due to the additional variable included in the ETM with respect to the STM (the volume fraction of blood plasma, vp ), the convergence of the optimization algorithm can be affected by the presence of local minima within the solution space. This means that the solution is dependent on the ini- tial seed. Although our model seems to overcome this issue on the simulated cases, the values obtained from the ETM could be used as an initial seed in those cases where the minimization convergence is severely affected by this issue. This new formulation is limited by the computational cost associated to the optimization algorithm. Although the method benefits from the computational efficiency of the FEM, the optimization time for the two simulated cases was 1.5 h and 10 h on average (fitting 360 time points on 339 nodes and 955 nodes, respectively). Despite being faster than the DP model (average of 72h to fit 484 voxels), the execu- tion time cannot be compared to the ETM, which took an average of 5s and 12s, respectively. Future works can be applied to migrate the code from Python to more efficient languages, such as C. One of the bottlenecks of this method is the Jacobian matrix computation, which executes opera- tions on large sparse matrix. This computation could be par - allelized to reduce the execution time.The described D-ETM is the first diffusion-corrected PK model to be implemented using the FE method. It pro- poses a new formulation for the diffusive term, based on the concept of effective diffusivity, that simplifies the com- putation of this term and avoids the inclusion of additional variables to the model. The semi-analytical method formu- lated to compute the Jacobian matrix opens the door for further gradient-based optimization methods for FE-based PK models. Although previous works [10] have developed a FE implementation of the extended Tofts model, their objective was to expose the effect of intra-voxel CA dif- fusion on PK analyses. To the best of our knowledge, this is the first FE-based optimization algorithm for the ETM. The results obtained with this model are promising, since it accurate retrieves the reference values, outperforming the ETM. Future works should test this model on real clinical or experimental data. Appendix A: Model parameters See Table 4. Note that although KTrans is measured in s−1 in the model, results are shown in min−1 to facilitate its comparison with the literature [8 , 15–18]. Appendix B: Numerical study of the diffusive term A small model (36 nodes and 120 timepoints) was created to compute its Jacobian matrix using the numerical approach based on the finite differences method, obtaining matrix JN , and the analytical method proposed in Sect. 2.3 (Eqs. 18, 19, 20), obtaining matrix JA . The numerical matrix JN was considered as ground truth. Given that most of values in the Jacobian matrix are zero or close to zero, we define a new error metric, ARDmax . This variable is defined as the quotient of the absolute difference between the Jacobian retrieved by each method and the maximum value of the numerical Jacobian. To improve the Table 4 List of model parameters and their unitsParameter Definition Units Ct Total CA concentration in the RVE mM/mm3 Ce CA concentration in the EES volume within the RVE mM/mm3 Cp CA concentration in blood plasma mM/mm3 KTrans Extravasation rate of CA from blood plasma to the EES 1/s ve Volume fraction of EES within the RVE 1 vp Volume fraction of plasma within the RVE 1 Deﬀ Effective diffusion coefficient mm2∕s D Diffusion coefficient in free medium mm2∕s 3863 Engineering with Computers (2022) 38:3849–3865 1 3 accuracy of the metric, we evaluate it for every nodal con- centration for each nodal variable and each timepoint (Eq. 22). E.g., for KTrans , we compute this metric for every nodal concentration for a certain timepoint ( ti ) for a certain nodal variable ( KTrans j): Figure 7 shows the error distribution of the analytical method with respect to the numerical one. The KTrans and vp Jacobian matrices show below 0.5% of error and 97% of the elements in the ve matrix show an error below 2%, what proves the accuracy of the analytical method.(22) ARDmaxti j=‖JAti j−JNti j‖ Max(JNti j)Appendix C: Benchmark problem: testing the influence of the axisymmetric geometry on convergence An additional simulation was generated to test the influence of the axisymmetric geometry on the convergence issues reported in Sect. 3.1. The geometry was modified to sup- press this symmetry, while the distributions of parameters were kept constant: KTrans = 0.2 min-1 in the rim and 0.05 min−1 in the core; a constant value of ve equal to 0.5 in the whole model; and vp=0.05 in the rim and 0.005 in the core. The initial seed was a random distribution of values similar to the one set in Sect. 3.1: between 0.4 min−1 and 0, for KTrans ; between 0.2 and 0.8 for ve ; and between 0 and 0.1 for vp . Ten different simulations were computed to evalu- ate the influence of the initial seed on the results. All of them converged to very similar values, validating the results obtained (Table 5). Despite that the D-ETM error the metrics in Table 5 show higher values than those obtained in Table 1, it still outperforms the ETM, especially on the KTrans map. The ETM, on the other hand, shows lower error in this case. This can be explained by the homogeneity in the parameters maps (Fig. 8). The reference maps in the benchmark prob- lem (Fig. 4) were heterogeneous, generating concentration gradients between adjacent elements within the vascular - ized region that produced diffusion fluxes. This diffusive process impacted on the parameters returned by the ETM, as explained in Sect. 3.1. In this case, the homogeneity on the parameters distributions avoids the influence of diffusion in the vascularized rim, what benefits the ETM. Fig. 7 Histogram of the metric ARDmax measured in the small model proposed. For KTrans and vp , the error is below 0.5%, while around 97% of the values for ve are below 2% Table 5 Comparison of error metrics between the D-ETM and the ETM for the non symmetric geometry of the benchmark problem The absolute error corresponds to the mean value of the model. KTrans absolute error is measured in min−1D-ETM ETM Absolute error Fraction of nodes error < threshold (%)Absolute error Fraction of nodes error < threshold (%) KTrans 1.14E-02 64.58 2.69E-02 45.03 ve 1.02E-01 56.14 1.73E+06 71.07 vp 8.45E-04 78.79 5.56E-03 77.71 3864 Engineering with Computers (2022) 38:3849–3865 1 3 Acknowledgements This work has received funding from the Euro- pean Union’s Horizon 2020 research and innovation programme under grant agreement No. 826494 (PRIMAGE). The work was also supported by the Ministry of Science, Education and Universities, Spain (FPU18/04541). Authors would like to acknowledge the Span- ish Ministry of Economy and Competitiveness through the project PID2020-113819RB-I00. The authors also would like to thank Silvia Hervas-Raluy, Francisco Merino-Casallo and Dr. Alberto Badias for their technical help and support. Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. Open Access This article is licensed under a Creative Commons Attri- bution 4.0 International License, which permits use, sharing, adapta- tion, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.