Front. Bioinform. Frontiers in Bioinformatics Front. Bioinform. 2673-7647 Frontiers Media S.A. 977228 10.3389/fbinf.2023.977228 Bioinformatics Original Research A proof-of-concept pipeline to guide evaluation of tumor tissue perfusion by dynamic contrast-agent imaging: Direct simulation and inverse tracer-kinetic procedures Vignon-Clementel et al. 10.3389/fbinf.2023.977228 Vignon-Clementel Irene E. 1 * Jagiella Nick 1 Dichamp Jules 1 Kowalski Jérôme 1 Lederle Wiltrud 2 Laue Hendrik 3 Kiessling Fabian 2 4 Sedlaczek Oliver 5 Drasdo Dirk 1 6 * 1 Inria, Palaiseau, France 2 Institute for Experimental Molecular Imaging (ExMI), University Clinic and Helmholtz Institute for Biomedical Engineering, RWTH Aachen University, Aachen, Germany 3 Fraunhofer MEVIS, Institute for Digital Medicine, Bremen, Germany 4 Fraunhofer MEVIS, Institute for Digital Medicine, Aachen, Germany 5 Department of NCT Radiology Uniklinikum/DKFZ Heidelberg, Heidelberg, Germany 6 IfADo - Leibniz Research Centre for Working Environment and Human Factors, Dortmund, Germany

Edited by: Jean-Karim Hériché, European Molecular Biology Laboratory Heidelberg, Germany

Reviewed by: Stefan Reinsberg, University of British Columbia, Canada

Mira Liu, University of Chicago Medicine, United States

 *Correspondence: Irene E. Vignon-Clementel, irene.vignon-clementel@inria.fr ; Dirk Drasdo, dirk.drasdo@inria.fr

This article was submitted to Computational BioImaging, a section of the journal Frontiers in Bioinformatics

 13 04 2023 2023 3 977228 24 06 2022 07 02 2023 Copyright © 2023 Vignon-Clementel, Jagiella, Dichamp, Kowalski, Lederle, Laue, Kiessling, Sedlaczek and Drasdo. 2023 Vignon-Clementel, Jagiella, Dichamp, Kowalski, Lederle, Laue, Kiessling, Sedlaczek and Drasdo

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Dynamic contrast-enhanced (DCE) perfusion imaging has shown great potential to non-invasively assess cancer development and its treatment by their characteristic tissue signatures. Different tracer kinetics models are being applied to estimate tissue and tumor perfusion parameters from DCE perfusion imaging. The goal of this work is to provide an in silico model-based pipeline to evaluate how these DCE imaging parameters may relate to the true tissue parameters. As histology data provides detailed microstructural but not functional parameters, this work can also help to better interpret such data. To this aim in silico vasculatures are constructed and the spread of contrast agent in the tissue is simulated. As a proof of principle we show the evaluation procedure of two tracer kinetic models from in silico contrast-agent perfusion data after a bolus injection. Representative microvascular arterial and venous trees are constructed in silico. Blood flow is computed in the different vessels. Contrast-agent input in the feeding artery, intra-vascular transport, intra-extravascular exchange and diffusion within the interstitial space are modeled. From this spatiotemporal model, intensity maps are computed leading to in silico dynamic perfusion images. Various tumor vascularizations (architecture and function) are studied and show spatiotemporal contrast imaging dynamics characteristic of in vivo tumor morphotypes. The Brix II also called 2CXM, and extended Tofts tracer-kinetics models common in DCE imaging are then applied to recover perfusion parameters that are compared with the ground truth parameters of the in silico spatiotemporal models. The results show that tumor features can be well identified for a certain permeability range. The simulation results in this work indicate that taking into account space explicitly to estimate perfusion parameters may lead to significant improvements in the perfusion interpretation of the current tracer-kinetics models.

 tumor perfusion in silico imaging tracer-kinetics models microcirculation dynamic contrast-enhanced imaging parameter estimation 864313 FKZ 0316042B European Research Council10.13039/501100000781 Bundesministerium für Bildung und Forschung10.13039/501100002347 1 Introduction

Certain histological features have been found to be characteristic for classification of tissue pathologies, in particular the grading of solid tumors and their evolution in absence of or in response to treatment (Noguchi et al., 2008; Hu et al., 2009; Li and Padhani, 2012; Essock-Burns et al., 2013), hence their knowledge can impact on clinical decisions. These features can be quantified by parameters such as the density and type of cells, the density, morphology and the spatial architecture formed by the vessels, which may also provide indirect information on leakiness (e.g. depending on fenestration of vessels etc.). However, histological information is not readily accessible unless biopsies are taken, which, besides being invasive, only provide local, and mainly morphological, information. For this reason, it is valuable to infer information on tissue microarchitecture from macroscopic imaging modalities, ideally complemented by functional signatures like the vessels’ flow and leakiness. Specific imaging modalities such as dynamic contrast-enhanced (DCE) perfusion imaging, or diffusion weighted imaging (DW)-MRI have been developed to gather information about such parameters down to the histology level non-invasively (Palmowski et al., 2008b; Lassau et al., 2011; O’Connor et al., 2011). As these modalities do not typically give direct access to the set of microscopic parameters, the signal measured with a given non-invasive imaging modality needs to be related to the information of interest at the histological level. This includes the spatial micro-architecture, as usually the signal corresponds to an average over many micro-structural elements, and the functional information as, for example, vessel leakiness, given that knowing the localisation of a blood vessel wall alone is insufficient to infer its permeability. Consequently, these non-invasive modalities, if correctly applied and interpreted, have the potential to provide histological and functional information, and, due to their non-invasiveness, can be applied during the development and treatment of the disease making them potentially very valuable for clinical decisions. For DW-MRI, the signal information can, under certain conditions, be easily related to the quantitative local tumour density (Yin et al., 2018) and even be probed by biopsies (Yin et al., 2021). For DCE imaging, numerous tracer kinetics (TK) models (e.g. Brix, Tofts, … ) have been constructed such that a small set of model parameters is believed to reflect the information of interest for the vasculature, and hence permit to retrieve this information from the imaging signal (DCE-MRI, DCE-CT, DCE-US).However, validation of these TK models remains a great challenge. Structural parameters can be compared to invasive histological data (Zwick et al., 2009). Functional parameters can be included for validation with in vitro phantom experiments (Driscoll et al., 2011; Gauthier et al., 2012). However to construct a representative vascularised tumor phantom at scale is challenging and difficult to adapt to different micro-architectures. Comparison to other blood flow imaging measurement modalities have validated flow (Lee et al., 2003), but only at the macro-scale. An alternative could be a “virtual” validation: by setting up a spatial temporal model that, even if it represents a simplification of the tissue micro-architecture, captures the important parameters and features of a tissue at the histological level, and then testing at least in silico, in how far TK-models can be expected to infer information on the parameters of interest, potentially proposing inference strategies that are expected to improve the inference of microscopic information from macroscopic signals.

We here establish a proof of concept for the entire path i.e., a workflow, starting with establishment of a direct model of blood vessel network and function, and simulating in how far two frequently used representatives of TK-models, e.g. the so called Brix and Tofts models, are able to infer information on tissue architecture and function from in silico DCE images (Figure 1).Our approach is not exhaustive but guided by Einstein’s principle that “everything should be made as simple as possible, but no simpler” to avoid unnecessary complexity. To limit the number of parameter sets, we focus on those parameter sets, that are capable of reproducing typical DCE-MRI signal patterns (Figure 2). These include the enhancement of a tumor’s border in case its microvascular density is elevated, while the core remains dark indicating necrosis (Figure 2A and movie TumorRimPerf.mov). Also an entire tumor can light up (in particular, if small) (Figure 2E) leading to characteristic fingerprints in the dynamic signal intensity (Figure 2F, f1, f2), which, for example, do not occur in muscle tissue (Figure 2F, f3). Even millimetre tumors are clinically observed (Figures 2C,D), which show a moderate elevated DCE signal (see also Supplementary Appendix SA8) and may display a necrotic core (Figure 2G).The workflow within this work connects a number of building blocks. For each of them, some aspects have already been addressed in published references, which for the sake of clarity we summarize only briefly in the introduction (more extensive information has been depicted in Supplementary Appendix SA2).

Benchmark concept and its different steps: in silico vascular network (A), in silico DCE-MRI generation (B), dynamic signal extraction (C), interpretation of the parameters with a tracer kinetics model (here a two-compartment model for extravascular contrast-agent (D)) by an inverse procedure (E), and comparison of obtained parameter maps (F) with the original ones. C A represents the input concentration, F blood flow, C P concentration in the vessels, V P their volume, K PS the exchange rate with the interstitial space, C I concentration in the latter, V I its volume.

Tissue perfusion characterization. Examples of tumors in human (upper row) and animal model (lower row). (A) DCE-MR-Image of Non-small cell lung cancer (NSCLC) tumor showing a characteristic enhanced ring after injection of a contrast agent (arrow: tumor of 4.1 cm (bar)). (B) Time course of the DCE-MRI signal intensity measured every 20s in a region of interest for an NSCLC tumor. (C) Prostate tumor (3.8 mm) in T2-weighted MRI. (D) DCE-MRI image of same tumor (circle). (E) DCE-derived parameter map of two Xenograft tumors (e1,e2) and muscle (e3) overlayed on T1-weighted MRI. (F) Corresponding signal intensity for tumors (f1,f2) and muscle (f3). (G) Histological cross-section of a few mm tumor grown in a Xenograft model showing cell nuclei (blue) surrounding a necrotic area (red).

The first contribution of this work is a flexible tool for vascular network creation and adaptation, control of features such as regional microvascular density and necrotic core, and computation for flow in this architecture (Figure 1A). The first step is flow modelling in blood vessel networks characterized by hallmarks of cancer vascularization, namely angiogenesis and necrosis (Hanahan and Folkman, 1996; Bergers and Benjamin, 2003; Vaupel, 2004). Stamatelos et al. (2014) were able to reconstruct the microvascular network of tumor Xenografts by a tremendous effort, but for demonstration of our workflow, the required network must be much greater than the tumor itself and its topology under control. Moreover, during the image reconstruction process, network connections may be missed out, which can be avoided if the network is constructed in silico. Representing hierarchical arteriovenous networks is however important to study tumor vascularization (Rieger and Welter, 2015). Hence by contrast to previous modelling work that mostly did not take it into account (Drasdo et al., 2010; Perfahl et al., 2011) or that rather focused on the development process (Welter et al., 2008; Welter and Rieger, 2010), we here concentrate on the generated tumor vascularization within a hierarchical network since it can be considered static at the time-scale of dynamic image acquisition. However, we only focus on the microvascular scale as these vessels mostly mediate the exchange of nutrients and contrast agents with the extravascular compartments and as they are representing the angiogenic sites.

The second contribution of this work is in the generation of in silico images from the coupled modeling of the contrast-agent in such a tumor vascular network and tissue, which here includes extravasation and diffusion in the interstitial space (Figure 1B). Indeed, given a computed blood flow, the transport of contrast agent or drugs inside these vessels and outside can then be computed. The spread of such compounds has been modelled either only intravascularly or considering a refined model of flow and drug concentration exchange, but not taking into account transport in the hierarchical networks feeding the tumor (Rieger and Welter, 2015). Closer to our work, Mescam et al. (2010) studied a multiscale model of blood flow and contrast-agent spread in liver, to better understand DCE-MRI images in the presence of a tumor. Flow is computed in realistic 3 d networks down to 40 μm connected with functional tissue units that do not resolve the microvascular (e.g. capillary) structure. By contrast, the model considered here simulates flow and transport in explicit vascular networks including capillaries, with transport computed inside the networks with emphasis on the essential exchange between the intra-and extravascular spaces. 2 and 3 d vascular networks as well as 1 d vessels are considered. The results show different functional examples and their influence on dynamic perfusion images.

Given a dynamic signal over time (Figure 1C), the next step is to infer some parameters from it, that reflect the underlying tissue (Figures 1E,F). The signal is interpreted at the image voxel level or in a whole region of interest (ROI), typically a region that is thought to be the tumor or a non-tumor zone, by a variety of TK models. None of the ones used in DCE imaging practice take space into account (Sourbron and Buckley, 2013). Semi-quantitative parameters (time to peak, maximum slope, etc.) can be easily extracted but without being specific for the different underlying tissue parameters. This work focuses on TK models that contain parameters surrogate for what is seen in typical tumors such as zones of higher microvascular density or significant leakage (Figure 2): the contrast agent concentration is often described by two interacting compartments, one for blood and one for the extravascular space (Tofts et al., 1999; Lee et al., 2003). It is assumed that the feeding concentration to the blood compartment is the same for all voxels, in some works with a delay-parameter (Lee et al., 2003; Koh et al., 2011; Sourbron and Buckley, 2013). Many TK models fit into the general theoretical framework proposed by Sourbron (2014).The third contribution of this work is in the assessment of the inverse procedure in 1, 2, and 3 d on synthetic DCE imaging data for two common TK models. Overall the goal is to validate the relationship between microscopic information and DCE imaging data in silico. Zwick et al. (2010) investigated how parameters of 2 TK models correlate to the tissue input parameters of an effective one-dimensional model of contrast-agent spread inside and out of the microvasculature i.e., without explicitly representing the vascular network and thus the tumor tissue structure in space. The same conceptual idea is pursued here, but starting from a more complex model, where space is explicitly taken into account both in the transport of the contrast agent in branching trees and connected tumor vascularization, and in the possible diffusion of the contrast agent outside the blood vessels, hence different from Zwick et al. (2010) generating spatial perfusion images and parameter maps (Figures 1B,D).

The structure of the paper (cf. Table of content before the Reference list) in Material and Methods and Results follows the different steps of the workflow (see Figure 1); namely functional microvascular network and tumor characteristics, generation of in silico dynamic perfusion images and estimation of the corresponding parameters from 2 TK models. This proof of concept of validating the relationship between microscopic data and DCE imaging data and the obtained results are finally discussed. The findings suggest that a more accurate parameter inference may be achieved if the distance of the feeding artery (if it is known) from the ROI is taken into account. The main components of each step of the workflow is explained in Section 2. This section also describes more precisely the novelties. Complementary technical information and results not necessary to understand the main conceptual line of the paper are given in the SI so that the paper is self-contained and the results can be reproduced (see also the Benchmark code section; some additional information can be found in ref (Jagiella, 2012)).

 2 Materials and methods

In the methods section, the spatial vascular architecture and perfusion model are described, followed by the contrast-agent intra- and extravascular transport modeling. Inference procedures are explained based on two common TK models similarly as in Brix et al. (2010).

2.1 <italic>In-silico</italic> tumor vascularization

In this subsection we briefly summarize qualitatively how the vasculature is generated with focus on the tumor vasculature (technical details, see Supplementary Appendix SB1). The construction of the vasculature requires and thereby implies the simulation of flow. The structure and remodeling of the microvascular networks closely follow previous works of Goedde and Kurz (2001) or Welter and Rieger (2010). But the purpose here is different: we consider the time scale of diagnosis, not of growth and remodeling. The vascularization is constructed to be directly reflective of different tumor architectures. Hence, novel local rules are proposed for the tumor region.

Vascular networks along with their flow properties (flow-rate, pressure, wall shear-stress) are constructed based on graphs placed on a regular lattice (in 1, 2, 3 d) by placing points into the center of each unit (a cube in 3d, square in 2 d, line segment in 1 d) and linking those with 3 d cylinders according to a certain set of rules specified in Supplementary Appendix SB1.1. Networks are first initialized from arterial and venous tree roots (Supplementary Appendix SB1.2), which become functional, interpenetrating trees connected by capillaries by shear-stress homogenization (Supplementary Appendix SB1.3). Vessel radii are recursively computed from capillaries up to the roots of these hierarchical networks, according to a power-law of coefficient α (equation (B.1)).

Inside growing tumors the micro-environment is different than in healthy tissue. In order to create a “tumor-like” vascularization, two regions - a tumor region and a necrotic core - are defined, where the rules and parameters for the homogenization algorithm are different (Supplementary Appendix SB1.3). Due to high proliferation, tumor cells consume much more nutrients than cells in the surrounding healthy tissue. Consequently, an increased production of VEGF in hypo-nourished parts of tumor stimulates vessel sprouting and increases microvascular density (MVD), and a delicate interplay of different growth factors and nutrient availability leads to different tumor vascularization scenarios (Drasdo et al., 2010; Welter and Rieger, 2010). During algorithmic vessel generation the parameter quantifying vessel sprouting (the “vessel sprouting probability” p tumor spr ) is increased while another parameter characterizing vessel collapse (the “vessel collapse probability” p tumor deg ) is decreased to reproduce a highly angiogenic tumor with PDGF (platelet-derived growth factor) mediated vessel maturation in the tumor model (see Table 1). In addition, the higher MVD is mimicked by connecting neighboring capillary ends with n parallel vessels, which modifies the relation of vessel radii at branching points and computation of flow along capillaries (Supplementary Appendix SB1.4).

Vascularization parameters specific to healthy, tumorous and necrotic tissue zones. τ denotes the wall shear stress, τ min and τ max minimal and maximal values, respectively (Supplementary Appendix SB1.3).

Normal tissue Tumor Necrotic core
Sprouting Probability p normal spr = 0.5 p tumor spr = 1 p necrotic spr = 0.1
Degeneration Probability p normal deg = τ max − τ τ max − τ min p tumor deg = p normal deg / 10 p necrotic deg = p normal deg ⋅ 10
Micro-Vessel Density n normal = 1 n tumor = 10 n necrotic = 10

The plasma volume fraction ϕ P , given by the ratio of volume occupied by vessels to total volume within a tissue region of interest, in tumorous tissue was found increased by one order of magnitude (ϕ P = 0.04 ± 0.01 in pectoral muscles; ϕ P = 0.2 ± 0.07 in carcinoma; see (Brix et al., 2004)). In the following, unless otherwise specified, the default value is n = 10 for inter-tip-connections inside tumoral regions (see Table 1).

Beyond a certain size, the central parts of the tumor are not nourished sufficiently anymore: a necrotic core appears. It was observed that this necrotic zone is not only lethal to tumor cells, but to endothelial cells and thus to small blood vessels as well (Kiessling et al., 2003; Palmowski et al., 2008a). Thus in this region, the vessel sprouting probability p necrotic spr is decreased and the vessel collapse probability p necrotic deg increased compared to healthy tissue (see Table 1).

 2.2 <italic>In-silico</italic> contrast-enhanced perfusion images

The vascular network structures are then used to simulate the transport of contrast-agents through vascularized tissues and create in silico perfusion images.

2.2.1 Contrast-agent modeling

Multi-phase equations model the macroscopic transport of the contrast agent in the tissue namely both inside the vascular network and the interstitial space. This model fits into the class of multiphase reaction-advection-diffusion models that for tracer kinetic models have been suggested to be generated by combining key process modules, namely, convection, diffusion, leakage, absorption, decay and sources (Sourbron, 2014). See Supplementary Appendix SB2.1 for its link to microscopic model and interphase exchanges.

The mass concentration of contrast-agent (CA) inside the plasma follows ∂ ∂ t ϕ P c P + ∇ ⋅ ϕ P c P v = − k P S c P − c I  in  Ω , c P ∂ Ω i n = c A t , where c A is the arterial (or more generally, vascular) input function (AIF), the blood plasma concentration in the arterial root nodes ∂Ω in . c P = n P /V P is the macroscopic (local averaged microscopic) plasma CA concentration, and c I = n I /V I is the macroscopic (local averaged microscopic) interstitial CA concentration, whereby n P , n I are the number of CA molecules in phases P and I, V P and V I the volumes of phases P and I in the tissue volume V = V P + V I . v is the macroscopic blood velocity, and k PS the volumetric macroscopic membrane exchange rate, both defined later. Parker et al. (2006) proposed an AIF functional form derived from patient data (see Supplementary Appendix SB.2.2, Supplementary Table SB4).

Assuming no advection but only diffusion in the interstitial space as a consequence of cells and extracellular matrix, and no-flux conditions at the tissue border, the CA mass concentration in the interstitial space follows ∂ ∂ t ϕ I c I = ∇ ⋅ ϕ I D I ∇ c I + k P S c P − c I  in  Ω , ∂ c I ∂ n ∂ Ω = 0 where ϕ I D I is the effective diffusion coefficient due to the fact that the molecules are not diffusing freely among the whole volume, but are limited by the cells and vessels borders. D I is taken as 103 μm 2/s (Lemke et al., 2009) while k PS will be defined later. For an in-depth derivation of the macroscopic equations from homogenization theory, see (Penta et al., 2015).

 2.2.2 Generation of parameter maps and <italic>in silico</italic> images

Eqs 1, 2 are then discretized and solved numerically for each voxel i (see Supplementary Appendix SB2.3). The flux f i,j from node i to node j relates to the velocity v i,j as f i , j = v i , j π r i , j 2 . The nodes are placed in the middle of each voxel (volume element on the regular lattice), i serves at the same time as a lattice identifier (Supplementary Appendix SB1.1). Parameter maps can then be defined for all voxels of an in silico image. Given the vascular graph G (defined in Supplementary Appendix SB1.1) gives for each voxel i its connected voxels j, the vascular volume fraction for voxel i of volume V i is given by ϕ P , i = V P , i V i = π 2 l 2 ∑ j∈ r i , j 2 . V P,i being the voxel vascular volume, and the extra-vascular volume fraction is ϕ I , i = 1 − ϕ P , i . The associated flow rate F i is the sum of all entering fluxes: F i = ∑ j ∈ connected neighbors f j , i | + . The volumetric exchange rate k PS,i in voxel i relates to the exchange rate K PS,i as k PS,i = K PS,i /V i in a voxel i. This exchange rate between the plasma and the interstitial space can be written as the product of the surface of the blood vessel walls S i,j , and the permeability coefficient, P i,j : K P S , i = ∑ j ∈ connected neighbors P i , j S i , j 2 = ∑ j ∈ connected neighbors P i , j π l r i , j By identification with the microscopic scale, P i,j represents D i , j m / H i , j m . Note, that half of the connecting vessel between two nodes i, j is associated with voxel i, the other with voxel j, which explains the division by a factor 2. In the rest of the paper, the permeability coefficient P i,j is a homogeneous constant noted P that can vary between healthy and tumor regions. Finally, the total concentration in a voxel i is then naturally defined as: c i t = ϕ P i c P i t + ϕ I i c I i t In case the ROI is a volume that contains several voxels, the in silico generated measurement is the total concentration averaged over all voxels of the ROI. The associated parameters are: ϕ P , R O I = ∑ i ∈ R O I ϕ P , i F ROI = ∑ i ∈ R O I ∑ j ∈ connected neighbors of i ∩ j ∉ R O I f j , i | + K P S , R O I = ∑ i ∈ R O I K P S , i Here it was used that all voxels are by construction of the same size, and K PS is an extensive quantity.Finally, the total concentration relates differently to the signal intensity depending on the image modality (Brix et al., 1991; Hoffmann et al., 1995; Brix et al., 2004). Without loss of generality, it is here directly taken as the surrogate for the in silico image signal over time in each voxel.

 2.2.3 Direct model parameters

The direct model parameters are not always known, and hence to choose representative values is a challenge. The network is constructed to reflect the increase from healthy to tumor tissue of plasma volume fraction (section 2.1). The healthy value is taken from muscle and can be changed. Tumor values are more variable (e.g. (Brix et al., 2004) reporting ranges of 0.1–0.4), and are extensively varied here (Supplementary Appendix SC3). A baseline permeability of p = 0.1 μm/s is used everywhere (in healthy and tumor zones). This value is inferred from (Brix et al., 2004; Brix et al., 2009) who found close K PS /(ϕ p *V) values for both muscle and carcinoma tissues, by estimating S/V and S/(ϕ p *V) from this paper’s networks and assuming a tissue density of 1g.mL −1. p = 0.1 μm/s was also reported for caco-2 (adenocarcinoma) and p = 0.35 μm/s MDCK (kidney) cells for paracetamol (Irvine et al., 1999), which has a diffusion coefficient of approximately 650 μm 2/s Ribeiro et al. (2012) close to what is chosen in this work. In Mescam et al. (2010); Zwick et al. (2010), K PS /V is varied from 0, 0.01, 0.1 and up to 1 min−1. Taking all this information into account motivated to vary P from 0, 0.01, 0.1–1 μm/s in most cases. These papers focused on gadolinium contrast agents, but permeability can be further varied.

 2.3 Estimation of functional parameters

DCE images give a value of the total CA concentration in every region of interest (ROI) for each acquisition time. The exact structure of the vascularization in each ROI and the interconnection between ROIs are completely unknown. In DCE imaging, simplified non-spatial models representing in each ROI separately the circulation of CA and its exchange with the rest of the tissue are fitted in an inverse procedure to recover perfusion and permeability information from the imaging dynamic curves.

2.3.1 Tracer kinetic models

In DCE imaging, a number of two-compartment models (Brix et al., 1991; Tofts and Kermode, 1991; Brix et al., 1999; Tofts et al., 1999; Sourbron, 2014) are the basis for quantification of regional blood flow, vessel permeability, and relative compartmental volumes within an ROI. In the following we will focus on two most common models proposed by Brix et al. (1999) and Tofts et al. (1999).

2.3.1.1 The Brix II model

This model, sometimes in the literature referred to as 2CXM (Sourbron and Buckley, 2013), takes into account in each ROI the blood (or plasma) and interstitial (extravascular) compartments (see Supplementary Figure SC13) of respective relative volumes ϕ P and ϕ I in ROI volume V, coupled by the intra-vascular-extra-vascular exchange rate K PS .

Each ROI blood compartment is fed directly from the AIF, with blood flow rate F going in and out of the compartment. No further transport or diffusion are taken into account. The corresponding equations are thus d ϕ P C P d t = F V C A − C P − K P S V C P − C I d ϕ I C I d t = K P S V C P − C I The total concentration of contrast agent in the tissue is thus assumed to be defined similarly to Equation 7. These equations are solved numerically with an implicit Euler scheme.

 2.3.1.2 The Extended Tofts model with and without delay

Tofts et al. (1999) proposed another two-compartment model (sometimes called Generalized Kinetic model), first neglecting the blood contribution assuming a negligible blood volume. d ϕ I C I d t = K trans C P − C I In the extended Tofts version, the blood pool contribution is added (Sourbron and Buckley, 2011). The total tissue concentration is then: C t = ϕ P C P t + K trans ∫ − ∞ t C P t ′ e − K trans / ϕ I ⋅ t − t ′ d t ′ . ϕ P and ϕ I being the plasma and interstitial compartment relative volumes respectively. The second term on the right hand side of Eq 12 denotes the solution of Eq 11, which emerges if the plasma volume is neglected (ϕ p = 0). C(t) is here explicitly given as a function of the CA plasma concentration C P (t) assumed to be the AIF, and of a transfer constant K trans × V whose meaning depends on the balance between the aforementioned flow F and permeability K PS (Tofts et al., 1999). As the AIF C A (t) might be measured in an area far away from the ROI, this assumption may not hold in most of the cases. Thus, a parameter t 0 can be added (Koh et al., 2011; Sourbron and Buckley, 2013), accounting for the time delay between the measured signal in the feeding artery and the vessels in the ROI. Eq. 12 becomes C t = ϕ P C A t − t 0 + K trans ∫ − ∞ t C A t ′ − t 0 e − K trans / ϕ I ⋅ t − t ′ d t ′ .

2.3.2 Parameter inference

The feeding concentrations C A (t) are assumed to be the same as the previously defined AIF c A (t). In practice, the dynamic signal for a given ROI is acquired with a certain temporal resolution Δt, for N t time steps of index n. This ROI-CA concentration at each time step c (n ⋅Δt) is compared to the TK solution C (n ⋅Δt) for this ROI with the L 2 norm S covering the entire acquisition time: S β = ∑ n = 1 N t c n ⋅ Δ t − C n ⋅ Δ t , β 2 Note that S(β) corresponds to the quadratic error. The same formula directly permits quantifying the goodness of the agreement between the ROI-CA and the TK solution of both Brix II and DE-Tofts model, as both have the same number of parameters (3, see below). To get the best fitting set of parameters β, S is minimized with the Levenberg-Marquardt algorithm with 100 iterations. For both tracer kinetics models, as a first approximation Eq. 4 is added as a constraint, since this is how the in silico dynamic signal is generated. For the Brix II model, β = {ϕ P , F, K PS } whereas for the Delay Extended Tofts Model (named thereafter DE-Tofts model) β = {ϕ P , K trans × V, t 0}.

An interesting information may be obtained by directly comparing the inferred parameters of the Brix II and the DE-Tofts model, which is straightforward for ϕ P and relatively straightforward for K PS = K trans × V, but there is no simple equivalence in the DE-Tofts model, that corresponds to the flow parameter F in Brix II. As the information of the volume flow rate F in the Brix II model has to be related to a parameter combination of the three Tofts model parameters, a first approach is by dimensional analysis, leading to the possible definition F = ϕ P V t 0 . However, as will be explained in the Results (Section 3.2.1), this definition needs to be modified. This modification for DE-Tofts gives a rational to do the same for the Brix II flow. When necessary to make a distinction between the flow F from the original map (Eq. 5) and the recovered one from the inference procedure, the latter is then labeled F recov .

 3 Results

First the results of the generation of in silico vascularizations and their validation against known architectural and functional properties are summarized. Then, the in silico dynamic contrast-enhanced images, the last step of the direct problem, are presented. The final part focuses on the results of the kinetics models assessment, i.e. the comparison of the parameters estimated with these models, against the original parameters from the in silico model. This leads to the proposal of a new interpretation of the estimated parameters to improve the perfusion assessment from the dynamic images.

Note that the 2 d models are 60 μm thick tissue slices, i.e. have a thickness of 60 μm ≪ extension in (x,y)-plane. The vascular graph itself lies on a 2 d plane with vessels having a certain cylinder volume, thus forming a connected network of a certain plasma volume as described in the Methods. The 3 d model is based on a vessel graph that extends in each of d = 3 dimensions. In all dimensions, each voxel has thus by default a size of (60 μm)3.

3.1 From vessel architecture to <italic>in silico</italic> dynamic contrast-enhanced images 3.1.1 Functional vascularization

Firstly, we verified that the vascular network generation algorithm (see Section 2.1) generates vessel networks that reproduce known morphological and functional features. Initially, the arterial and venous trees stay well separated from each other (Supplementary Figure SC11). Large parts of the vascular networks have low to almost no flow, and, consequently, are not functional. The construction algorithm is iterative; with the number of iterations, inter-tree-connections increase and with them the overall flow, resulting in interpenetrating hierarchical arterio-venous networks with a homogeneous wall shear stress in both the 2d and 3 d models (see Supplementary Appendix SC1).

Importantly, the model results actually compare well with measured data on wall shear stress, mean vessel segment flow, pressure-velocity distribution and pressure, in terms of order of magnitudes and changes with vessel radius (Figure 3, right column).

Comparison with experimental measurements. Left: simulation results for three different arterial-venous trees with pairs of exponent α (denoted by A for the arterial, V for the venous tree), with the same applied inlet and outlet pressures: in black the pair (A,V)=(3,2), in red the pair (2.7,2.4) and in blue the even pair (2.7,2.7). For shear stress, velocity and pressure, arteries are shown on the negative diameter side, while veins are on the positive diameter side. For flow their trees are plotted together to demonstrate that flow in arteries and veins follow the same power laws as the corresponding diameter (d α ), and can be identified by their exponent. Right: literature data. The different rows show wss (from Fung (2013)), flow (from Huo and Kassab (2007)), velocity (from Pries and Secomb (2008)) and pressure (from Fronek and Zweifach (1974)), versus diameter.

The precise relations depend on the exponents α art and α ven of the power-law (equation B.1). In several statistics of the obtained micro-circulation properties, we observe differences between vascular trees (e.g. pressure, see Figure 3) and asymmetries between venous and arterial vessels (e.g. wall shear stress and blood velocity, see Figure 3). Three different in silico vascularizations in 2d, with varying exponents α art and α ven are considered. The higher the exponent α, the higher the flow for a given applied pressure drop (Figure 3, second row) indicating a drop of flow resistance. This results from two competing effects that can also explain the differences of decay between flow rate, velocity and wall shear stress with respect to diameter when increasing the exponent (see Supplementary Appendix SC2). The best fit of the asymmetric properties between model and data leads to choose for the rest of the article different default values for arteries, α art = 3 and for veins, α ven = 2.7 (Figure 3).

 3.1.2 Vascularized tumor cases

After having generated normal vascularizations, we now study vascularized tumor cases for various MVDs reflected in the parameter value of n defined in equation B.5) (Supplementary Figure SC12). Above a certain MVD volume fraction ϕ P becomes locally sufficiently different from its values in neighboring regions (second row in Supplementary Figure SC12). The higher volume fraction leads to a higher flow, also well detectable already at smaller n (last row in Supplementary Figure SC12). In the necrotic core (Supplementary Figure SC12 last column), the model induces the destruction of the small vessels in it, but the larger vessels remain inside it. In the tumor periphery, there is an increased ϕ P . The latter corresponds to a characteristic observation in patient tumors (Figure 2A). This structure of larger vessels remaining inside tumors with low vessel density in the center and high density at the periphery is typically observed in multi-modality imaging in Xenograft experiments (Kiessling et al., 2004b,a). The later also reports dilated vessels feeding the tumor as emerges in the model here as well.

 3.1.3 <italic>In-silico</italic> dynamic contrast-enhanced images

Based on microcirculation reflective of the normal and tumor regions, in silico DCE images (quantifying c(t) as defined in Eq. 7), illustrate the kinetics specific to each component (Figure 4; and movies for further varied parameters in Supplementary Appendix SC3). The reference case leads to a quite homogeneous vascular volume fraction and hence capillary-tissue exchange rate K PS . The larger vessels are visible due to their larger values. The corresponding in silico overall image intensities follow in their time course the kinetics of the AIF, with the peak of the first pass, then a large decrease, then a second peak followed by a slower decay. Compared to the average image intensity, the larger vessels are highlighted, and the zones of capillary-tissue exchange appear as a haze. By contrast, as soon as MVD is locally increased (n = 30) in the tumor region, this increases the local capillary-exchange rate. On top of the previously described dynamics, the tumor outer-rim always appears highlighted. Modelling the tumor by an increased permeability (p = 100 instead of 0.1 μm/s) only also leads to a local capillary-tissue exchange rate increase. However, the in silico DCE-images are different (Figure 4 bottom row): the tumor is highlighted much later, with a diffuse glow around the entire tumor zone. The explanation is that in case of elevated MVD the simulated signal intensity reflects directly the MVD and hence, follows the kinetics of tumor perfusion, dictated by the AIF, while in case of an elevated permeability, the marker leaks into the extra-vascular space and returns as soon as the vessel marker concentration drops below the extra-vascular marker concentration, which generates the delay.

Increased MVD and/or permeability. Increased MVD in the tumor region (first row) by comparison to the reference simulation (second row), with permeability p = 0.1 μm/s for these two simulations. Increased permeability in the tumor region (p = 100 μm/s, last row) by contrast to the reference simulation, with normal MVD for these two simulations. Plasma volume fraction (first column, ϕ p ), exchange rate K PS (second column, in μm 3/s) and in silico DCE images at four different times (remaining columns). K PS in μm 3/s. Voxels of (60 μm)3, domain:100 × 100 voxels. Movies: see list in Supplementary Appendix SC3.

These results correspond to well-known tumor enhancement differences. For example, hepatocellular carcinomas, neuroendocrine tumors, renal cell carcinomas, and squamous cell carcinomas have strong arterial feeding and high MVD and thus, are clinically characterized by their early enhancement after contrast agent administration in CT and MRI scans. The other extreme are low to moderately differentiated gliomas where the disruption of the blood brain barrier leads to enhanced vascular permeability. These tumors are best seen in late enhancement phases after contrast agent injection, when the contrast agent had sufficient time to slowly distribute in the extravascular tumor space. These special vascular features are of high importance for diagnosing and differentiating many tumors and have strongly influenced the clinical guidelines on the performance of CT and MRI examinations.

Hence, the different components of what defines a vascularised tumor lead to different dynamical features of the DCE-images. It is thus natural to extract quantitative differences in parameters based on the dynamics differences seen in DCE images. This is the focus of the next section.

 3.2 Inference results

The problem of inferring model parameters from a dynamic intensity in a region of interest is called an inverse problem, or inverse procedure. This procedure requires to have a measurement (here the concentration over time) and to select a model to interpret this measurement. In-vivo the exact mathematical model that represents this dynamics is not known, and in practice simplified models are selected. Here we generate in silico images with a known model, and perform the inverse problem with simplified models (Brix II and DE-Tofts models) to better understand their relevance.

3.2.1 Intravascular agents

As a reference case we study intravascular agents in one, two and three dimensions (= d) (Figure 5). The advantage of performing in silico simulations is the opportunity of such a strategy starting with a great simplification, and then step-wise studying more complex cases gives the opportunity of gaining a step-wise understanding.

Inference for zero permeability (A) 1 d case geometry (green voxel of (60 μm)3) and flow parameter along the 100-voxel domain (data, recovered and rescaled). (B) 2 d case of 100 × 100 voxels with the data (left column) and the Brix II recovered maps (right column) (DE-Tofts: see (C) 15) for (60 μm)3 voxels. The feeding artery is in the center of the image (c1): vascularization topology (from high (red) to low (blue) pressure) and four different regions of interest (ROI size = 709 voxels, ≈ 2 mm diameter). I contains the tumor necrotic zone (red circle) inside the tumor rim (green circle). Feeding artery in the center of the image (c2): the fit (continuous curves) of the measured contrast-agent concentrations over time (circles) for the four zones by the Brix II model. (D) Parameter inference from 3 d data (6mm × 3mm × 3 mm). The feeding artery is on the left of the images (central slices in z-direction) (d1) ϕ p (1rst row) and flow (second row) for a (60 μm)3 resolution; (d2) ϕ p (1rst row) and flow (second row) for a 300 × 300 × 3000 μm 3 resolution. Columns: data, Brix II (A–D) F is in μm 3/s.

 3.2.1.1 The one dimensional case

In a first step we study the inverse procedure based on Brix II and Tofts models in the simplemost case, which is a one-dimensional straight piece of tissue, as for this case the concentration over time can be analytically calculated. Despite this case looks extremely simplified, it may be compared to the case where the ROI is small compared to the distance between the feeding point and the ROI, in which case the Euclidean distance between feeding point and ROI may be of the order of the distance that the blood has to travel between them. The tissue piece is assumed to be a rectangular cuboid of length of 6 mm and cross-sectional area of L*L = 60 μm*60μm, containing a straight cylindrical capillary inside of constant radius r = 4 μm (thus ϕ P ≈ 0.014) along the long axis (see Supplementary Figure SC13). An intravascular agent enters at the inlet (at x = 0) with the known concentration C A (t) taken equal to the AIF. In that case, the plasma concentration c P is obtained from Eq. (7) by dropping the vessel leakage term to c P x , t = C A t − x / v = C A t − x π r 2 / F = C A t − x ϕ P L 2 / F , representing “travelling wave” i.e., a concentration profile that performs a rigid movement with a certain velocity v = F/(ϕ P L 2) (F is the volume flow rate, ϕ P the plasma volume fraction).

The common case is when the ROI does not contain the feeding point of the contrast-agent (Figure 5A1, Supplementary Figure SC13C) (the case where the ROI contains the feeding point is considered in Supplementary Appendix SC4, Supplementary Figure SC13A, B). Consider as a measurement a fixed volume V i = L 2*l i = (60 μm)3 but which center is at varying distance x from the feeding point (green volume in Figure 5A1). Considering conservation of mass, equation C.1 becomes d c ̄ d t = d ( ϕ P c ̄ P ) d t = F V i ( C A ( t − ( x − l i / 2 ) / v ) − c P ( x + l i / 2 , t ) ) = F V i ( C A ( t − ( x − l i / 2 ) / v ) − C A ( t − ( x + l i / 2 ) / v ) ) . In this case the ratio of volume and flow, or equivalently l i /v, are small relative to the characteristic duration of the signal represented by the time width of the CA concentration peak (Figure 1), l i /v < time-width of CA peak, or equivalently F/(ϕ P V i ) > 1/time-width of CA peak. The fact of replacing the outlet concentration at position x i + l i /2 by the average concentration in the inverse Brix II model (Eq. (9)) is a good approximation. The issue is the inlet flux: with increasing x, taking C A (t) instead of C A (t − (x − l i /2)/v) is becoming decreasingly accurate (see Supplementary Appendix SC4.1). Thus Eq. (9) under-predicts more and more the flow, as the estimated flow F recov ≈ F/N i , where N i = x/l i (Figure 5A2) i.e., F recov x is a constant that depends on the chosen size of the ROI, l i . This is validated by the fact that rescaling the recovered flow F recov by multiplying it by N i , one recovers the exact flow (Figure 5A2, Supplementary Figure SC13C).

In a next step, we consider the DE-Tofts model for the same configuration (Supplementary Figure SC13D). In this model, we have apriori defined the estimated flow as F recov = ϕ P V/t 0. In our simple case, Eq. 13 written for an intravascular agent reduces to ϕ P C A (t − t 0). On the other hand, the analytical solution (Eq. 15) at x identifies t 0 = xϕ P L 2/F. So for an ROI small enough, in the sense described above of F/(ϕ P V) > 1/time-width of CA peak, the concentration is sufficiently homogeneous in the ROI and l ≈ N i *l i . Hence when taking this ROI concentration for the parameter inference, one recovers t 0, and thus F recov = ϕ P V/t 0. But as in fact, t 0 ≈ N i ϕ P V/F, the recovered flow F recov ≈ F/N i . Thus, here too it needs to be rescaled by the voxel-distance to the feeding point to recover the original flow F.

Sensitivity to noise: As described in Supplementary Appendix SC4.2, a measurement error or random noise in the tissue measurement affects the inference of the flow more than the inference of the plasma volume fraction. The latter is always estimated with an accuracy of ∼±20% even for noise up to 20%.

 3.2.1.2 Two dimensional case

In a next step the case of intravascular agent propagation through a vasculature which contains a vascularised tumor is considered (Figure 5B). The in silico contrast-agent images are generated as in section 3.1 (displayed in Supplementary Figures SC11, 12), with blood and the contrast-agent entering through an artery in the middle of the tissue, the contrast-agent feeding temporal profile following the typical AIF of equation Supplementary Appendix SB2.2. The corresponding plasma volume fraction (Figure 5B1) and flow (Figure 5B2) range from 0 to 0.14 and 0–106  μm 3/s respectively, while P = 0 μm/s. The measured concentration in each voxel, computed as the solution to Eq. (7), is the input to the inverse procedure.

Without assuming that the agent is intravascular, the estimated parameters are thus the plasma volume fraction, the plasma-tissue exchange coefficient and the flow. For Brix II, the estimated plasma volume fraction map is visually identical to the initial one (Figure 5B1 right vs. Figure 5B1 left): the main vasculature structure, and the vascularized tumor zone are well detected. Finally, the estimated flow map (Figure 5B2 right), shows the microvasculature structure and the vascularised tumor zone. Due to the rescaling of the flow with the Euclidean distance of the voxel to that feeding point, the flow maps are close to the initial data (Figure 5B2 left). Without the rescaling, the flow going away from the central feeding point is rapidly underestimated, similar to the 1d case (Supplementary Figure SC15). The DE-Tofts model gives very similar results to the Brix II model (Supplementary Figures SC15 right). While in regions where the plasma volume fraction is non-zero, the exchange coefficient is negligibly small as expected, in regions of zero plasma volume (and very low flow), as for example in the necrotic region of the tumor, the exchange coefficient deviates from zero (Supplementary Figures SC15 third row).

 3.2.1.3 ROIs of multiple voxels

Next the signal considered for the inverse problem is for a region of interest (ROI) larger than one voxel, as this is often the case in practice (Figure 5C). First, a 2 d vascularised tumor case is modelled with the transport of an intravascular contrast agent (Figure 5C1). The measurements from which parameters are estimated are the spatial average concentrations in four different ROIs, three in the normal zone (II, III, IV) and one containing a tumor (I). Each region contains 709 voxels, thus are all of diameter around 2 mm (typically a few in-plane voxels on DCE-MRI images). The concentrations for the four ROIs are similar to the AIF, although with a peak scaling from 0.46 to 1 mM and some time-dispersion due to ROI averaging (see zoom picture in Figure 5C2). For both Brix II and DE-Tofts, the measurement curves are well reproduced after parameter estimation, for all four regions, with small deviations between the two first passes as seen in the zoom pictures (Figure 5C2, Supplementary Figures SC16A3). The inverse procedure leads to estimated parameters, one for each ROI, that are compared to the original ones (as defined over each ROI according to Eqs 1, 2) in Table 2 The plasma volume fractions range from 0.08 for zone IV to 0.2 for zone I, overestimated on average for both models by 27%. The rescaled flows range from 5 107  μm 3/s for zone IV to 108  μm 3/s for zone I, with an average error for both models of 13%. In the inference procedure, the plasma-tissue exchange parameter was not set a priori to zero, and was recovered low for both models (at least four order of magnitude lower than the flow). Despite these errors, the four zones are ranked in the same order as the original data in terms of plasma volume fraction and flow, for both models: the region with the highest flow and plasma volume fraction is the one that contains the tumor.

Data (direct problem) and estimated parameters for Brix II and DE-Tofts for each ROI I, II, III and IV of Figure 5. K PS and F are in μm 3/s. ROI size: 709 voxels, ≈ 2 m m diameter and ≈ 0.15 m m 3 volume V. Relative errors ((data-recovered)/data) are reported for ϕ P and F. For K PS , since data is 0, no relative error can be reported.

Zone Data Brix II DE-tofts Error brix II Error DE-Tofts
ϕ P K PS F ϕ P K PS F ϕ P K PS F ϕ P F ϕ P F
I 1.67E-01 0 1.26E+08 2.05E-01 2.71E+03 1.08E+08 2.00E-01 9.30E+03 1.12E+08 0.23 0.14 0.20 0.11
II 7.61E-02 0 6.24E+07 9.78E-02 1.90E+02 8.27E+07 9.70E-02 1.20E+03 8.16E+07 0.29 0.33 0.27 0.31
III 6.63E-02 0 5.24E+07 8.50E-02 9.48E+02 5.15E+07 8.35E-02 3.04E+03 5.32E+07 0.28 0.02 0.26 0.02
IV 6.10E-02 0 4.73E+07 8.12E-02 3.80E+02 5.09E+07 7.99E-02 1.96E+03 5.09E+07 0.33 0.08 0.31 0.08
3.2.1.4 The three dimensional case

Finally, a 3 d example with an intravascular agent is presented, where the vascularized tumor tissue is built in the same manner as for the previous 2 d case but with now a 3 d network graph (Figure 5D). In this example, the vascular tumor rim is a thick sphere: 2 d central cuts of the resulting flow and plasma volume fractions show the main vessels, the tumor rim and necrotic core. When the inverse procedure is performed for each voxel (as in Supplementary Appendix SC5), only Brix II gives reasonable results, i.e. the recovered parameter maps show the same structures as the original ones (Figure 5D1, right vs. Supplementary Figures SC16B1 right). Despite the Euclidean distance rescaling of the flow, one can see the influence of the distance to the feeding artery in the quality of the results. For the DE-Tofts model, this influence gets rapidly worse (Supplementary Figures SC161, right). In such conditions, considering ROIs of the size of a typical DCE-MRI voxel (0.3 mm*0.3 mm*3 mm) leads to a major blur of the parameters (Figure 5D2, Supplementary Figures SC16B2). In the original data, the larger vessels and tumor zone cannot be numerically distinguished in the flow, but only the tumor zone in the plasma volume fraction can. The same patterns are seen in the Brix II recovered parameters. As expected by the results at the voxel level, the DE-Tofts model barely recovers the tumor zone in plasma voume fraction map (Supplementary Figures SC16B2, right).

In summary, the simulation results show that plasma volume fraction, exchange rate and the flow can be reasonably well estimated, the latter by scaling the flow rate with the Euclidean distance between feeding point and ROI. This tendency holds true in d = 1, two for the Brix II and DE-Tofts models. In 3 d the Brix II model performs better than the DE-Tofts model. Generally, the quality of parameter inference decreases with d.

 3.2.2 Extravascular agents

Next, extravascular cases are considered. Direct simulations for four different levels of permeabilities P are carried out in one and quasi-3d cases, to test the influence of permeability level on the recovery procedure. For the quasi-3D cases, the initial maps are inherently composed of varying vessel density, flow and exchange parameters. Each time the inference procedure is tested for every voxel with the 2 TK models.

3.2.2.1 The one dimensional case

Different from the purely intravascular situation (Figure 5A) now the exchange rate can be non-zero. The default value is taken as P ∼ 0.1 m/s and in the simulations a range of permeability values is considered (see Section 4.3 on the choice of these values). In the direct simulations, the permeability is varied by several orders of magnitude, resulting in the CA tissue concentration deviating more and more from the AIF, both close and far from the inlet (Supplementary Figures SC17): the first peak gets lost in the signal that continues to rise much later than the AIF, and significantly higher than the first peak contrarily to the AIF. The different models can in general well capture these shapes. The plasma volume fractions and exchange rates (Supplementary Figures SC18 bottom and middle respectively) are well estimated with the different models, but this becomes increasingly more difficult with increasing permeability and distance to the inlet, eventually failing for P = 1 μm/s. The apparent flow parameters estimated by the different models differ from the input flow. Flows can be recovered as in the previous section by rescaling with the length of the path from the feeding artery (Supplementary Figures SC18). But this is getting more inaccurate as permeability increases. In fact, the models become unable to capture the signal shape for slow apparent flows typically recovered in voxels far from feeding artery (Supplementary Figures SC17 third row) for P = 1 μm/s. Because of boundary conditions at the entrance and exit of the tube, an artefact can occur for these voxels (Supplementary Figures SC17).

 3.2.2.2 The multi-dimensional case

Next the 2 d case is considered as in Figure 5B, but this time for an extravascular agent. Note that even though the vasculature graph is embedded in 2 d (plane), the simulation is 3 d in the sense that diffusion in the extravascular space occurs everywhere in the domain, like in a 3 d slab of tissue (2 d plane that includes the vasculature over a 6 mm*6 mm square and a small height of 60 μm). We refer to this case as ‘quasi-3d’. Different permeabilities are modelled keeping the other parameters unchanged. Thus the flow and plasma volume fraction parameter maps are the same for p = 0.1 μm/s (Figure 6) and P = 1 μm/s (Figure 7) as for the intra-vascular agent described above (Figure 5, Supplementary Figures SC15). The resulting plasma-tissue exchange coefficient by contrast increases with an increasing permeability, making apparent the necrotic zone (very low K PS ), the larger blood vessels and the tumor rim (very high K PS ). The measured concentration c(t) in each voxel, solution to Eqs 1, and, is then the input to the inverse procedure. We do not show p = 0.01 μm/s, because the results are as close to P = 0 μm/s as in the 1 d case. For p = 0.1 μm/s (Figure 6 two last columns), the inferred plasma volume fraction is well captured by both models. The rescaled flow and plasma-tissue exchange coefficient are also well recovered, except in zones where there are no blood vessels. For P = 1 μm/s (Figure 7 two last columns), the situation is markedly different: BrixII is the only one to recover the structures of the flow maps but not the other two parameters, whereas with DE-Tofts, the map structures of plasma volume fraction and K PS (where there are blood vessels) are captured but the flow structures are very degraded.

Inferred TK model parameters from contrast-agent time course in quasi-3d. Columns: left = data, center = Brix II, right = DE-Tofts, vessel permeability is p = 0.1 μm/s. First row: vascular density ϕ P , second row: flow, third row: rescaled flow, fourth row: exchange rate K PS . Voxel size: (60 μm)3, domain: 100 × 100 × 1 voxels. F and K PS in μm 3/s. Same scale bar for all ϕ P , F, K PS .

Inferred TK model parameters from contrast-agent time course in quasi-3d. Columns: left = data, center = Brix II, right = DE-Tofts, vessel permeability is P = 1 μm/s. First row: vascular density ϕ P , second row: flow, third row: rescaled flow, fourth row: exchange rate K PS . Voxel size: (60 μm)3, domain: 100 × 100 × 1 voxels. F and K PS in μm 3/s. Same scale bar for all ϕ P , F, K PS .

 4 Discussion

The main contribution of the paper is the proof of concept presented in Figure 1, which is based on spatial in silico perfusion images. We constructed in silico vasculatures to solve the direct problem of contrast agent perfusion (intra-vascular transport, intra-extra-vascular exchange and diffusion within interstitial space) for evaluation of parameter estimations from different tracer kinetics models used in dynamic contrast-enhanced (DCE) perfusion imaging. The transport in and out of the 2 d or 3 d vascular networks is modeled by multiphase equations that use as an input a 1 d model of blood flow inside the vessels, the latter limiting the computational cost. Here, as a proof of principle we showed the evaluation procedure of Brix II and DE-Tofts with in silico perfusion data of a contrast agent after a bolus injection. For the inverse procedure, six ground truth 1D cases (varying permeability and noise levels), four ground truth 2D/quasi-3D cases (varying permeability levels), and one ground truth 3D case were carried out as direct simulations. For 2D/3D cases, vessel volume fraction, flow and exchange parameters were inherently varying locally. Different inverse procedures were tested each time - i.e. varying ROI location and size, and TK models. For each 2D simulation, the inverse problem was performed for the 104 voxels, and in 3D, the inverse problem was performed for 2.5 105 voxels and for the coarser case for 200 ROIs.

4.1 Discussion on the generated architecture and DCE-images

The direct model parameters are not always known (see section 2.2.3), and hence to assess if they are representative of real tissues is a challenge. In healthy vasculature (e.g. Figure 6), K PS = 200 μm 3/s, and thus K PS /(ϕ p *V) = 1.4 min−1, which is close to the value of 0.2 min−1 reported in (Brix et al., 2004). In the tumor region in Figure 6 where parameter inference is successful for p = 0.1 μm/s, K PS /V = 0.3 min−1, which is close to the values reported by (Brix et al., 2009) and used as baseline in (Zwick et al., 2010) K PS /V = 0.2 ml min−1. g −1. Regarding flows, the obtained healthy and tumor flows F/V ≈ 15–150 min−1 on a voxel level and F/V ≈ 20–40 min−1 for the ROI cases (Table 2) are higher than the values reported in (Brix et al., 2004; Brix et al., 2010) but the increase by a factor of 2–10 from healthy to tumor values is reasonable. Moreover, a number of healthy perfusion parameters compared well with literature data over several vessel diameters and tree architectures (Section 3.1.1).

For the tumor vasculature and DCE-images, qualitative comparisons were made in Sections 3.1.2, 3.1.3 to further assess the realism of simulations. Besides, the change of shape of the concentration over time, when increasing permeability (e.g. Supplementary Figures SC17) is similar to that reported in Figure 2B for a lung tumor or also in (Brix et al., 2004; Zwick et al., 2010).

The model can be readily adapted if more precise parameter information is available, concerning each component of the pipeline, including the vascular network generation and their related functional values (if e.g. more detailed experimental studies simultaneously providing information on micro-vascular network architecture and function become available).

 4.2 Discussion on the parameter inference assessment

Regarding the interpretation of the perfusion parameters, overall parameter inference through TK models successfully detected the tumor area, with increased vascular density, flow and exchange rate (leakage). Studying the influence of space, we thus propose the following parameter inference procedure:

• directly interpret the plasma volume fraction,

• where the former is non-negligible, consider the plasma-tissue exchange-rate recovered value,

• rescale the flow-rate by the distance of the voxel to the contrast-agent arterial input function location.

Similarly to (Zwick et al., 2010), we found that the vessel volume fraction was a very robustly recovered parameter in all cases: with or without noise, in multiple dimensions, varying permeability. For the perfusion parameter (flow) however, we showed how results from DCE-imaging parameter estimation need to be interpreted with care. This is especially true the further the voxel is from the feeding artery: in this case, the flow is more underestimated compared to its true value.

We proposed to rescale the value with the distance to the feeding artery. Actually, the real distance traveled by the contrast-agent should be taken into account, so by rescaling with the true graph distance (if accessible) instead of by the Euclidean distance to the feeding artery the parameter inference would be expected to further improve. A simple estimate of this distance may be obtained by computing the average time of transport through the identified network and the blood flow, yielding a flow-related distance. The distributed parameter models based on the 1 d cylindrical capillary-tissue system are attempts to take this spatial component into account but they have not gained practical use (Koh et al. (2011); Sourbron and Buckley (2013)), probably because the overall network leads to more complex spatial influences (see discussion on dispersion below). A TK model that take only diffusion fluxes into account but with a known diffusion parameter, has been shown to improve the inverse procedure compared to the standard Tofts model for a direct model that does not included advection but focuses on diffusion (Pellerin et al., 2007).

For permeability or the plasma-tissue exchange rate K PS , the recovered value was reasonable if the intrinsic permeability P is not too high compared to F/S i.e., if the dimensionless number P/(F/S) = (PS)/F = K PS /F was not too high (see for example Supplementary Figures SC17). The value should be discarded though when the vessel volume fraction or perfusion are close to zero, because in this case the signal intensity is not sensitive to this parameter. For the DE-Tofts model, this is consistent with the interpretation of K trans as being K PS only in high flow regime (Tofts et al., 1999).

The 2 d simulations (Figures 6, 7, quasi-3d in that they have been performed in a 2 d slice of thickness 60 μm) suggest a possible strategy by combining both BrixII and DE-Tofts models to infer the parameter set composed of plasma volume, plasma-tissue exchange rate and flow (even if unusual for DE-Tofts). If both give similar results for the plasma volume, the permeability is likely to be not too high, and both models may be used to estimate all three parameters. In our simulations, very similar results were obtained for zero up to moderate plasma-tissue exchange rates. If both give different values, this may indicate a high permeability. In that case, BrixII gives reasonable estimates of the flow and plasma-tissue exchange rate, while DE-Tofts better captures the plasma volume fraction. However, a complex 3 d network architecture fed by a distant feeding vessel (as in Figure 5) surprisingly hindered a reliable parameter estimation in particular by the DE-Tofts model even for zero vessel permeability, indicating further need of extensive studies for 3 d vasculatures.

This work also highlights the dispersion effects introduced by heterogeneous vessel networks. Although interesting (Bassingthwaighte and Goresky, 1984), this is another challenge for parameter estimation: the current tracer kinetics models do not contain dispersion per say, and thus the recovered parameters are affected by it. Nevertheless, in 2d, the ranking of the different zones (more or less vascularized and with functional flow) could be distinguished. More work is needed in this direction. In particular, one could use the direct problem to model dispersion in a discretized vascular network in a similar approach as Koch et al. (2020).

In fact, as the size of the ROI increases, dispersion and other effects e.g. due to several feeding vessels entering and leaving the ROI, might affect the overall ROI signal. As a consequence, the recovered ROI parameters will be effective, and may or may not relate to the microscopic parameter as in Eq. 8. The preliminary results on the 3 d case indicate that varying ROI sizes or image resolution impacts the structures that can be identified. To study how ROI parameters scale with ROI size will allow to compare their values with those reported in vivo at typical image resolution. In that sense, the proposed framework can be used to study the effects of the perfusion image resolution on understanding the underlying microscopic structure and function.

 4.3 Current limitations

Several simplifications have been made in the current framework that could be improved as future work. First, the vascularization is built on a regular lattice and thus the vessel diameter is limited by the lattice constant. This constant is 60 μm, reflecting a reasonable orthogonal minimal capillary distance of the order of the diffusion length of oxygen in tissue and typical values determined in tissue micro-architectures but the precise value would depend on the particular tissue under consideration (Pries and Secomb, 2011). To consider larger vessels and thus larger trees, generation of the vasculature in lattice-free space could be favorable even if computationally more expensive. Here the choice was made to have tumors embedded into an otherwise normal hierarchical vasculature, and as a consequence these tumors are small. Although tumors of similar sizes (2–4 mm) are found histologically and on DCE-MRI data (Figures 2C,D,G), studying larger tumors with the same proof of concept would be a valuable future complement to the current work.

Furthermore, the extravascular compartment does not distinguish between intracellular and extracellular spaces, to reduce the number of estimated parameters. As long as the cells do not take up the CA as well, this is not changing the qualitative kinetics meaning only the precise concentration values may change. I.e. it does not affect the conclusions of this work. In all cases the same workflow described in this work can be pursued to study the precise parameter values of the different inference models.

Finally, additional physical effects could be implemented in the direct problem such as tumor filtration as done by Penta et al. (2015), without changing the workflow.

 4.4 Strengths of approach

The model is flexible since one can vary architecture and function of the vascularization (control of features such as regional microvascular density, necrotic core, etc.) and see their effects in synthetic images. This framework allows to test hypotheses supported by that another recent work points into similar directions (Nuha Abdul et al., 2021). Different other TK models can also be tested within the same framework. Moreover, such model permits to study the ambiguity caused by different tumor vasculatures (highly permeable versus high microvessel density) leading to similar perfusion images (bright tumor area, dark necrotic core) and their characteristic differences (temporal occurrence as e.g. delayed clearance in case of high permeability). In that sense, it brings function to histological data. While histological slices (see an example in Figure 2G) provide direct tissue information, they only give vessel density information or, after challenging image analysis, architecture, but not functional information. Moreover, often part of the vessels found in pathological specimen are not perfused. In-silico modeling may thus provide an interesting complement to interpret such data.

Different injection procedures (bolus injection versus gradual administration) could be easily studied along with other implications for diagnosis (temporal resolution and duration of acquisition).

 5 Conclusion

To conclude, this work presents the construction of a pipeline to investigate the link between DCE images and the underlying tissue microstructure and function, and demonstrates its possibilities on a variety of cases. This pipeline is based on a spatial model of vascularization and contrast-agent transport in and out of the vasculature. The direct model parameters for healthy and tumor zones are based on literature values and varied for sensitivity analysis of the in silico images to model inputs. Overall the parameter inference from the TK models typically applied to characterize DCE images, successfully detects the tumor features compared to the surrounding tissue. The results of the parameter inference indicate that the inference for both Brix II and delay-adjusted extended Tofts models is robust for low permeability and little dispersion due to network heterogeneity. The plasma volume fraction and the vascular-tissue exchange rate where the former one is non-negligible, are directly accessible while the recovered flow parameter needs some reinterpretation. Rescaling the flow by the ROI distance to the feeding artery significantly improves its estimation, the Euclidean distance being a work-around up to a certain point to estimate the time the injected CA volume needs to travel to the ROI. The results indicate the importance of taking explicitly into account space in inverse models. This study lays the foundation for better studying the effect of coarse-graining (or ROI sizes) in such dynamic imaging modalities and also for better interpreting histological data since they only provide structural but not functional parameters.

 Data availability statement

The original contributions in the study are publicly available. This data can be found here: https://gitlab.inria.fr/simbiotx/dyna-imaging-mod.

 Author contributions

In decreasing order of contribution: Conceptualization: DD, IV-C, OS, FK Methodology: IV-C, NJ, HL, DD Software: NJ, JK, JD Validation (e.g. patient and xenograft examples): OS, FK, HL, WN Visualization: NJ, JD Writing/original draft: IV-C Writing/review and editing: DD, IV-C, JD, FK, HL, OS Project admin and funding: DD, IV-C.

We acknowledge the funding sources of the German Ministry of Education and Research through LungSys (FKZ 0316042B) as well as LiSyM-Cancer (FKZ 031L0257D) and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 864313).

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbinf.2023.977228/full#supplementary-material

References Bassingthwaighte J. Goresky C. (1984). “Modeling in the analysis of solute and water exchange in the microvasculature,” in Handbook of physiology. The cardiovascular system. Microcirculation (New York: American Physiology Society), 4, 549–626. Bassingthwaighte J. (2021). Mmid4. Available at: http://www.physiome.org . Bassingthwaighte J. Wang C. Chan I. (1989). Blood-tissue exchange via transport and transformation by capillary endothelial cells. Circulation Res. 65, 997–1020. 10.1161/01.res.65.4.997 Bergers G. Benjamin L. E. (2003). Tumorigenesis and the angiogenic switch. Nat. Rev. Cancer 3, 401–410. 10.1038/nrc1093 Boissier N. Drasdo D. Vignon-Clementel I. E. (2021). Simulation of a detoxifying organ function: Focus on hemodynamics modeling and convection-reaction numerical simulation in microcirculatory networks. Int. J. Numer. Methods Biomed. Eng. 37, e3422. 10.1002/cnm.3422 Brix G. Bahner M. L. Hoffmann U. Horvath A. Schreiber W. (1999). Regional blood flow, capillary permeability, and compartmental volumes: Measurement with dynamic CT—initial experience. Radiology 210, 269–276. 10.1148/radiology.210.1.r99ja46269 Brix G. Kiessling F. Lucht R. Darai S. Wasser K. Delorme S. (2004). Microcirculation and microvasculature in breast tumors: Pharmacokinetic analysis of dynamic MR image series. Magnetic Reson. Med. 52, 420–429. 10.1002/mrm.20161 Brix G. Salehi Ravesh M. Griebel J. (2017). Two-compartment modeling of tissue microcirculation revisited. Med. Phys. 44, 1809–1822. 10.1002/mp.12196 Brix G. Semmler W. Port R. Schad L. R. Layer G. Lorenz W. J. (1991). Pharmacokinetic parameters in CNS Gd-DTPA enhanced MR imaging. J. Comput. Assisted Tomogr. 15, 621–628. 10.1097/00004728-199107000-00018 Brix G. Zwick S. Griebel J. Fink C. Kiessling F. (2010). Estimation of tissue perfusion by dynamic contrast-enhanced imaging: Simulation-based evaluation of the steepest slope method. Eur. Radiol. 20, 2166–2175. 10.1007/s00330-010-1787-6 Brix G. Zwick S. Kiessling F. Griebel J. (2009). Pharmacokinetic analysis of tissue microcirculation using nested models: Multimodel inference and parameter identifiability. Med. Phys. 36, 2923–2933. 10.1118/1.3147145 Cai Y. Wu J. Li Z. Long Q. (2016). Mathematical modelling of a brain tumour initiation and early development: A coupled model of glioblastoma growth, pre-existing vessel Co-option, angiogenesis and blood perfusion. PLoS One 11, e0150296. 10.1371/journal.pone.0150296 Cattaneo L. Zunino P. (2014). A computational model of drug delivery through microcirculation to compare different tumor treatments. Int. J. Numer. methods Biomed. Eng. 30, 1347–1371. 10.1002/cnm.2661 Courant R. Friedrichs K. Lewy H. (1928). Über die partiellen Differenzengleichungen der mathematischen Physik. Math. Ann. 100, 32–74. 10.1007/bf01448839 Drasdo D. Jagiella N. Ramis-Conde I. Vignon-Clementel I. E. Weens W. (2010). Modeling steps from a begnin tumor to an invasive cancer: Examples of instrinsically multiscale problems. United Kingdom: Taylor & Francis, 379–416. Driscoll B. Keller H. Coolens C. (2011). Development of a dynamic flow imaging phantom for dynamic contrast-enhanced CT. Med. Phys. 38, 4866–4880. 10.1118/1.3615058 Essock-Burns E. Phillips J. J. Molinaro A. M. Lupo J. M. Cha S. Chang S. M. (2013). Comparison of DSC-MRI post-processing techniques in predicting microvascular histopathology in patients newly diagnosed with GBM. J. Magnetic Reson. Imaging 38, 388–400. 10.1002/jmri.23982 Fick A. (1855a). Ueber diffusion. Ann. Phys. Chem. 170, 59–86. 10.1002/andp.18551700105 Fick A. (1855b). V. On liquid diffusion. Philos. Mag. Ser. 4 10, 30–39. Forsgren M. F. Leinhard O. D. Dahlström N. Cedersund G. Lundberg P. (2014). Physiologically realistic and validated mathematical liver model revels hepatobiliary transfer rates for Gd-EOB-DTPA using human DCE-MRI data. PLoS One 9, e95700. 10.1371/journal.pone.0095700 Fronek K. Zweifach B. W. (1974). Pre- and postcapillary resistances in cat mesentery. Microvasc. Res. 7, 351–361. 10.1016/0026-2862(74)90022-3 Fung Y. c. (2013). Biomechanics: Circulation. Berlin, Germany: Springer Science & Business Media. Fung Y. C. (1990). Biomechanics: Motion, flow, stress, and growth. New York: Springer-Verlag. Gabrani-Juma H. Clarkin O. Pourmoghaddas A. Driscoll B. Wells G. deKemp R. (2016). Validation of a multimodality flow phantom and its application for assessment of dynamic spect and pet technologies. IEEE Trans. Med. Imaging 36 (1), 132–141. 10.1109/TMI.2016.2599779 Gauthier M. Pitre-Champagnat S. Tabarout F. Leguerney I. Polrot M. Lassau N. (2012). Impact of the arterial input function on microvascularization parameter measurements using dynamic contrast-enhanced ultrasonography. World J. radiology 4, 291. 10.4329/wjr.v4.i7.291 Goedde R. Kurz H. (2001). Structural and biophysical simulation of angiogenesis and vascular remodeling. Dev. Dyn. 220, 387–401. 10.1002/dvdy.1118 Hanahan D. Folkman J. (1996). Patterns and emerging mechanisms of the angiogenic switch during tumorigenesis. Cell. 86, 353–364. 10.1016/s0092-8674(00)80108-7 Hoffmann U. Brix G. Knopp M. V. He T. Lorenz W. J. (1995). Pharmacokinetic mapping of the breast: A new method for dynamic MR mammography. Magnetic Reson. Med. 33, 506–514. 10.1002/mrm.1910330408 Hu L. Baxter L. Smith K. Feuerstein B. Karis J. Eschbacher J. (2009). Relative cerebral blood volume values to differentiate high-grade glioma recurrence from posttreatment radiation effect: Direct correlation between image-guided tissue histopathology and localized dynamic susceptibility-weighted contrast-enhanced perfusion MR imaging measurements. Am. J. Neuroradiol. 30, 552–558. 10.3174/ajnr.a1377 Huo Y. Kassab G. S. (2007). A hybrid one-dimensional/Womersley model of pulsatile blood flow in the entire coronary arterial tree. Am. J. Physiol. Heart Circ. Physiol. 292, H2623–H2633. 10.1152/ajpheart.00987.2006 Ingrisch M. Sourbron S. (2013). Tracer-kinetic modeling of dynamic contrast-enhanced MRI and CT: A primer. J. Pharmacokinet. Pharmacodyn. 40, 281–300. 10.1007/s10928-013-9315-3 Irvine J. D. Takahashi L. Lockhart K. Cheong J. Tolan J. W. Selick H. E. (1999). MDCK (Madin-Darby canine kidney) cells: A tool for membrane permeability screening. J. Pharm. Sci. 88, 28–33. https://onlinelibrary.wiley.com/pdf/10.1021/js9803205. Jagiella N. (2012). Parameterization of lattice-based tumor models from data. Ph.D. thesis (Paris, France: Université Pierre et Marie Curie-Paris VI). Kiessling F. Farhan N. Lichy M. P. Vosseler S. Heilmann M. Krix M. (2004a). Dynamic contrast-enhanced magnetic resonance imaging rapidly indicates vessel regression in human squamous cell carcinomas grown in nude mice caused by VEGF receptor 2 blockade with DC101. Neoplasia 6, 213–223. 10.1593/neo.03394 Kiessling F. Greschus S. Lichy M. P. Bock M. Fink C. Vosseler S. (2004b). Volumetric computed tomography (VCT): A new technology for noninvasive, high-resolution monitoring of tumor angiogenesis. Nat. Med. 10, 1133–1138. 10.1038/nm1101 Kiessling F. Heilmann M. Vosseler S. Lichy M. Krix M. Fink C. (2003). Dynamic T1-weighted monitoring of vascularization in human carcinoma heterotransplants by magnetic resonance imaging. Int. J. cancer 104, 113–120. 10.1002/ijc.10913 Koch T. Flemisch B. Helmig R. Wiest R. Obrist D. (2020). A multiscale subvoxel perfusion model to estimate diffusive capillary wall conductivity in multiple sclerosis lesions from perfusion MRI data. Int. J. Numer. Methods Biomed. Eng. 36, e3298. 10.1002/cnm.3298 Koh T. S. Bisdas S. Koh D. M. Thng C. H. (2011). Fundamentals of tracer kinetics for dynamic contrast-enhanced MRI. J. Magnetic Reson. Imaging 34, 1262–1276. 10.1002/jmri.22795 Kurz H. Sandau K. (1997). Modelling of blood vessel development - bifurcation pattern and hemodynamics, optimality and allometry. Comm. Theor. Biol. 4, 261–291. Lassau N. Chami L. Chebil M. Benatsou B. Bidault S. Girard E. (2011). Dynamic contrast-enhanced ultrasonography (DCE-US) and anti-angiogenic treatments. Discov. Med. 11, 18–24. Lee D. S. Rieger H. Bartha K. (2006). Flow correlated percolation during vascular remodeling in growing tumors. Phys. Rev. Lett. 96, 058104. 10.1103/physrevlett.96.058104 Lee T. Y. Purdie T. G. Stewart E. (2003). CT imaging of angiogenesis. Q. J. Nucl. Med. 47, 171–187. Leedale J. A. Kyffin J. A. Harding A. L. Colley H. E. Murdoch C. Sharma P. (2020). Multiscale modelling of drug transport and metabolism in liver spheroids. Interface Focus 10, 20190041. https://royalsocietypublishing.org/pdf/10.1098/rsfs.2019.0041. Lemke A. Laun F. B. Klauß M. Re J. T. Simon D. Delorme S. (2009). Differentiation of pancreas carcinoma from healthy pancreatic tissue using multiple b-values: Comparison of apparent diffusion coefficient and intravoxel incoherent motion derived parameters. Investig. Radiol. 44, 1–7. Li S. P. Padhani A. R. (2012). Tumor response assessments with diffusion and perfusion MRI. J. Magnetic Reson. Imaging 35, 745–763. 10.1002/jmri.22838 Luypaert R. Sourbron S. de Mey J. (2011). Validity of perfusion parameters obtained using the modified Tofts model: A simulation study. Magn. Reson Med. 65, 1491–1497. 10.1002/mrm.22728 Macklin P. McDougall S. Anderson A. R. Chaplain M. A. Cristini V. Lowengrub J. (2009). Multiscale modelling and nonlinear simulation of vascular tumour growth. J. Math. Biol. 58, 765–798. 10.1007/s00285-008-0216-9 McDougall S. R. Anderson A. R. Chaplain M. A. (2006). Mathematical modelling of dynamic adaptive tumour-induced angiogenesis: Clinical implications and therapeutic targeting strategies. J. Theor. Biol. 241, 564–589. 10.1016/j.jtbi.2005.12.022 Mescam M. Kretowski M. Bezy-Wendling J. (2010). Multiscale model of liver DCE-MRI towards a better understanding of tumor complexity. IEEE Trans. Med. Imaging 29, 699–707. 10.1109/tmi.2009.2031435 Noguchi T. Yoshiura T. Hiwatashi A. Togao O. Yamashita K. Nagao E. (2008). Perfusion imaging of brain tumors using arterial spin-labeling: Correlation with histopathologic vascular density. Am. J. Neuroradiol. 29, 688–693. 10.3174/ajnr.a0903 Nuha Abdul U. N. Hanson E. A. Sourbron S. Munthe-Haas A. Z. (2021). A computational fluid dynamics framework to generate digital reference objects for perfusion imaging objects for perfusion imaging. ISMRM & SMRT (ANNUAL MEETING, ABSTRACT No 0497), 4. O’Connor J. Tofts P. Miles K. Parkes L. Thompson G. Jackson A. (2011). Dynamic contrast-enhanced imaging techniques: CT and MRI. Br. J. radiology 84, S112–S120. 10.1259/bjr/55166688 Palmowski M. Huppert J. Hauff P. Reinhardt M. Schreiner K. Socher M. A. (2008a). Vessel fractions in tumor xenografts depicted by flow- or contrast-sensitive three-dimensional high-frequency Doppler ultrasound respond differently to antiangiogenic treatment. Cancer Res. 68, 7042–7049. 10.1158/0008-5472.can-08-0285 Palmowski M. Huppert J. Ladewig G. Hauff P. Reinhardt M. Mueller M. M. (2008b). Molecular profiling of angiogenesis with targeted ultrasound imaging: Early assessment of antiangiogenic therapy effects. Mol. Cancer Ther. 7, 101–109. 10.1158/1535-7163.mct-07-0409 Parker G. J. M. Roberts C. Macdonald A. Buonaccorsi G. A. Cheung S. Buckley D. L. (2006). Experimentally-derived functional form for a population-averaged high-temporal-resolution arterial input function for dynamic contrast-enhanced MRI. Magnetic Reson. Med. 56, 993–1000. 10.1002/mrm.21066 Pellerin M. Yankeelov T. E. Lepage M. (2007). Incorporating contrast agent diffusion into the analysis of DCE-MRI data. Magnetic Reson. Med. 58, 1124–1134. 10.1002/mrm.21400 Penta R. Ambrosi D. Quarteroni A. (2015). Multiscale homogenization for fluid and drug transport in vascularized malignant tissues. Math. Models Methods Appl. Sci. 25, 79–108. 10.1142/S0218202515500037 Perfahl H. Byrne H. M. Chen T. Estrella V. Alarcón T. Lapin A. (2011). Multiscale modelling of vascular tumour growth in 3D: The roles of domain size and boundary conditions. PloS one 6, e14790. 10.1371/journal.pone.0014790 Pinder G. F. Gray W. G. (2008). Essentials of multiphase flow in porous media. Wiley-Interscience. Pries A. R. Secomb T. W. (2008). “Blood flow in microvascular networks,” in Handbook of physiology: Microcirculation. Editors Tuma R. F. Duran W. N. Ley K. Second edition (San Diego: Academic Press), 3–36. Pries A. R. Secomb T. W. (2011). Blood flow in microvascular networks. Atlanta, Georgia, United States: American Cancer Society, 3–36. chapter 1. Pries A. R. Secomb T. W. Gessner T. Sperandio M. B. Gross J. F. Gaehtgens P. (1994). Resistance to blood flow in microvessels in vivo . Circ. Res. 75, 904–915. 10.1161/01.res.75.5.904 W. K. Purves D. Sadava G. H. Orians H. C. Heller (Editors) (2003). Life: The science of Biology. 7th edition (W.H.Freeman and Co Ltd.), 954. Ribeiro A. C. Barros M. C. Veríssimo L. M. Santos C. I. Cabral A. M. Gaspar G. D. (2012). Diffusion coefficients of paracetamol in aqueous solutions. J. Chem. Thermodyn. 54, 97–99. 10.1016/j.jct.2012.03.014 Rieger H. Welter M. (2015). Integrative models of vascular remodeling during tumor growth. Wiley Interdiscip. Rev. Syst. Biol. Med. 7, 113–129. 10.1002/wsbm.1295 Secomb T. W. Alberding J. P. Hsu R. Dewhirst M. W. Pries A. R. (2013). Angiogenesis: An adaptive dynamic biological patterning problem. PLoS Comput. Biol. 9, e1002983. 10.1371/journal.pcbi.1002983 Sourbron S. (2014). A tracer-kinetic field theory for medical imaging. IEEE Trans. Med. imaging 33, 935–946. 10.1109/tmi.2014.2300450 Sourbron S. P. Buckley D. L. (2013). Classic models for dynamic contrast-enhanced MRI. NMR Biomed. 26, 1004–1027. 10.1002/nbm.2940 Sourbron S. P. Buckley D. L. (2011). On the scope and interpretation of the Tofts models for DCE-MRI. Magn. Reson Med. 66, 735–745. 10.1002/mrm.22861 Sourbron S. P. Buckley D. L. (2012). Tracer kinetic modelling in MRI: Estimating perfusion and capillary permeability. Phys. Med. Biol. 57, R1–R33. 10.1088/0031-9155/57/2/r1 Stamatelos S. K. Kim E. Pathak A. P. Popel A. S. (2014). A bioimage informatics based reconstruction of breast tumor microvasculature with computational blood flow predictions. Microvasc. Res. 91, 8–21. 10.1016/j.mvr.2013.12.003 Tofts P. S. Brix G. Buckley D. L. Evelhoch J. L. Henderson E. Knopp M. V. (1999). Estimating kinetic parameters from dynamic contrast-enhanced t1-weighted MRI of a diffusable tracer: Standardized quantities and symbols. J. Magnetic Reson. Imaging 10, 223–232. 10.1002/(sici)1522-2586(199909)10:3<223:aid-jmri2>3.0.co;2-s Tofts P. S. Kermode A. G. (1991). Measurement of the blood-brain barrier permeability and leakage space using dynamic MR imaging. 1. Fundamental concepts. Magnetic Reson. Med. 17, 357–367. 10.1002/mrm.1910170208 Vaupel P. (2004). Tumor microenvironmental physiology and its implications for radiation oncology. Semin. Radiat. Oncol. 14, 198–206. 10.1016/j.semradonc.2004.04.008 Welter M. Bartha K. Rieger H. (2008). Emergent vascular network inhomogeneities and resulting blood flow patterns in a growing tumor. J. Theor. Biol. 250, 257–280. 10.1016/j.jtbi.2007.09.031 Welter M. Rieger H. (2013). Interstitial fluid flow and drug delivery in vascularized tumors: A computational model. PLoS One 8, e70395. 10.1371/journal.pone.0070395 Welter M. Rieger H. (2010). Physical determinants of vascular network remodeling during tumor growth. Eur. Phys. J. E Soft Matter 33, 149–163. 10.1140/epje/i2010-10611-6 Wu M. Frieboes H. B. Chaplain M. A. J. McDougall S. R. Cristini V. Lowengrub J. S. (2014). The effect of interstitial pressure on therapeutic agent transport: Coupling with the tumor blood and lymphatic vascular systems. J. Theor. Biol. 355, 194–207. 10.1016/j.jtbi.2014.04.012 Yang X. Knopp M. V. (2011). Quantifying tumor vascular heterogeneity with dynamic contrast-enhanced magnetic resonance imaging: A review. J. Biomed. Biotechnol. 2011, 732848. 10.1155/2011/732848 Yin Y. Breuhahn K. Kauczor H. U. Sedlaczek O. Vignon-Clementel I. E. Drasdo D. (2021). Diffusion-weighted mri-guided needle biopsies permit quantitative tumor heterogeneity assessment and cell load estimation. arXiv: 2103.00714. Available at: https://arxiv.org/abs/2103.00714 (preprint). Yin Y. Sedlaczek O. Müller B. Warth A. González-Vallinas M. Lahrmann B. (2018). Tumor cell load and heterogeneity estimation from diffusion-weighted MRI calibrated with histological data: An example from lung cancer. IEEE Trans. Med. Imaging 37, 35–46. 10.1109/tmi.2017.2698525 Zwick S. Brix G. Tofts P. S. Strecker R. Kopp-Schneider A. Laue H. (2010). Simulation-based comparison of two approaches frequently used for dynamic contrast-enhanced MRI. Eur. Radiol. 20, 432–442. 10.1007/s00330-009-1556-6 Zwick S. Strecker R. Kiselev V. Gall P. Huppert J. Palmowski M. (2009). Assessment of vascular remodeling under antiangiogenic therapy using DCE-MRI and vessel size imaging. J. Magn. Reson Imaging 29, 1125–1133. 10.1002/jmri.21710 Nomenclature V

Graph vertex set

 E

Graph edge set

 G = (V, E)

Graph

 l i,j

Length of edge (i, j)

 r i,j

Radius of edge (i, j)

 η(r i,j )

Apparent viscosity of edge (i, j)

 G i , j = π 8 ⋅ r i , j 4 η ( r i , j ) ⋅ l i , j

Hydraulic conductance of edge (i, j)

 Δp i,j = p i − p j

Pressure difference on edge (i, j)

 f i,j = G i,j ⋅Δp i,j

Flow rate on edge (i, j)

 τ i , j = 4 η ( r i , j ) f i , j π r i , j 3

Shear stress on edge (i, j)

 ϕ P , i = V P , i V i

Plasma volume fraction in voxel of volume V i

ϕ I,i = 1 − ϕ P,i

Extra-vascular volume fraction in voxel of volume V i

F

Flow rate in voxel of volume V i (volume flow rate)

 c m

Microscopic concentration

 D m

Diffusion coefficient

 H m

Membrane thickness

 S m

Membrane surface

 J = −D m ∇c m

Diffusion flux

 v m

Microscopic velocity field of the fluid containing the contrast agent

 c P

Plasma macroscopic concentration

 k PS

Volumetric macroscopic membrane exchange rate

 k P S i

Volumetric macroscopic membrane exchange rate in a voxel i

K P S i

Exchange rate in a voxel i

P i , j = D i , j m H i , j m

Permeability of edge (i, j)

 P

Permeability constant over a region for all its vascular edges

 v

Macroscopic velocity field

 c A , C A

Arterial input function concentration

 c I

Extra-vascular macroscopic concentration

 D I

Diffusion coefficient in extravascular domain

 c i = ϕ Pi c Pi + ϕ Ii c Ii

Total concentration in voxel i

C P

Brix II/Extended Tofts plasma concentration

 C I

Brix II/Extended Tofts extravascular concentration

 K trans

Transfer constant in extended Tofts model

 t 0

Tofts delay

 β

Set of parameters to fit

 S(β)

Quadratic fitness function to minimize

 F recov

Estimated flow rate from compartment models