Perfusion estimation using synthetic MRI-based measurements and a porous media flow model

The measurement of perfusion and filtration of blood in biological tissue give rise to important clinical parameters used in diagnosis, follow-up, and therapy. In this paper, we address techniques for perfusion analysis using processed contrast agent concentration data from dynamic MRI acquisitions. A new methodology for analysis is evaluated and verified using synthetic data generated on a tissue geometry.


Introduction
Accurate estimates for blood perfusion are of high importance in clinical applications.Perfusion imaging is used in treatment planning for stroke [1], but also for clinical applications related to the staging of cancer [2].
Medical image acquisition techniques like computerized tomography (CT), magnetic resonance imaging (MRI), or positron emission tomography (PET), can all be applied in a dynamic setting where the evolving distribution of an injected contrast agent is gathered together as a temporal sequence of images.Quantitative tissue characterization (e.g., blood perfusion) from such data is traditionally performed locally by applying tracer-kinetic methodology [3] to a single region of interest (ROI) or a voxel at a time.While these approaches have advantages in terms of computational cost, they are generally considered to have certain systematic weaknesses and the resulting estimates will typically lack robustness in terms of accuracy and reproducibility, see e.g., [4][5][6].
In the current work, we clearly distinguish between flow and perfusion and understand perfusion to be the delivery of blood to the tissue, i.e. the capillary system.This is necessary to account for both larger-scale transport of blood in an organ and micro-scale utilization of the blood.The traditional models [3] do not have a spatial component.When applied to spatialtemporal data one thereby implicitly assume well-defined boundaries and known contrast inflow distribution given as an arterial input function (AIF) for any ROI.This is a valid assumption if whole organs or large ROIs are considered but do not necessarily hold when smaller ROIs, such as voxels are investigated.The observed contrast agent concentration for a voxel or a small group of voxels will typically reflect an aggregation of flow in transit in the arterial and/or venous system as well as flow relevant to the local tissue perfusion or extravascular leakage.Local temporal deconvolution also neglects the potentially useful spatial structure of the images.Advanced fluid models and numerical simulators developed for subsurface flow can be used to predict the blood flow in organs [7].Models for flow in reservoir rocks have a large number of poorly known parameters.This is also the case when modeling blood flow in organs, as exact tissue properties are impossible to know, and are also subject to large variations between different individuals [8].Ensemble-based data assimilation is a popular approach to improve models for fluid flow, based on available measurements (see e.g., [9]).Iterative ensemble smoothers (see e.g., [10,11]) have become the preferred approach for assisted history matching of petroleum reservoir properties.Smoothers assimilate all data in one single update step, and complicated and time-consuming restarts of the simulator are avoided.However, a problem with ensemble-based methods is that with a limited ensemble size, spurious correlations between observations and model parameters are inevitable.This problem increases with the size of the parameter space and the number of observations and leads to unrealistic updates of model parameters, and underestimation of the posterior error covariance, see [12].As a remedy for this, localization may be introduced in order to limit updates of parameters that we know are unlikely to influence a given observation.In our study we suggest applying correlation-based localization, [13], where a threshold for the spurious correlations between the model parameters and the concentration measurements are estimated using the ensemble at the current iteration.
In this paper, we investigate the potential of using fluid flow models in combination with an ensemble smoother for improved estimation of blood perfusion from dynamic MRI.Unknown parameters in a porous media flow model (based on Darcy's equation) are updated, using a time sequence of spatially distributed contrast concentration data.Synthetic data are generated using a hybrid-scale model where flow in the largest vessels is simulated using the Hagen-Poiseuille equation for laminar incompressible flow and the tissue flow is modeled using Darcy's equation.This hybrid-scale model is more complex compared to the pure porous media model, and the goal is to reconstruct the blood flow and perfusion despite the gap in complexity between the two models.The synthetic data have known, ground truth, perfusion maps.The methodology is tested on a domain representing a stretched-out frog tongue, with known vascular structure.A comparison is made with the classical maximum slope method and deconvolution based on singular value decomposition.

Methods
We introduce two different modeling approaches.One to be used in parameter estimation and a more complex and computationally demanding to be used to generate synthetic data with known ground truth.For the purpose of data assimilation, we consider a slightly modified version of the traditional dual-porosity, dual permeability model (see e.g., [14]).This partitioning of the geometry into multiple compartments allows an explicit characterization of the local tissue perfusion as the local transfer flux between the compartments.Thus enabling a clear distinction between fluid contributing to local perfusion and fluid in transit through the voxel.
Second, we formulate a hybrid model where a Darcy-type porous media model is combined with two network models representing arterial and venous flow respectively.This model will be used in lieu of real observations to provide tracer-concentration data for the assimilation exercise as well as ground truth perfusion maps.The hybrid model is developed and designed to run on a fine grid and the equations are solved using an implicit solution strategy.This allows for larger time steps which again is important to reduce the computational time.

Porous media model used for parameter estimation
We model a region of live tissue as a three-dimensional spatial domain O partitioned into four compartments.Three of the compartments constitute the vasculature and are all represented as porous media flow models.The fourth compartment represents the extravascular tissue, and will for now be considered a passive tissue matrix confining our pore spaces.
Let the subscript i 2 {a, v} indicate the arterial or venous compartments respectively.As detailed in [7], we consider two porous flow models each governed by the Darcy constitutive relation u i = −μ −1 K i (x)rp i , the incompressibility condition r � u i = q i , and conservation of contrast agent Here, u is Darcy velocity, μ is blood viscosity, K is the permeability tensor, p is pressure, q is the source (if positive) or sink (if negative) term, ϕ is porosity, and c is concentration.The departing (q < 0) fluid carries the overall concentration of the compartment, i.e., c * a ¼ c a , while in general the entering (q > 0) fluid carries a different concentration, i.e., c * v 6 ¼ c v .The concentration carried by the entering fluid c * v equals c * a shifted by a capillary transit time.This procedure is explained in detail below.The arterial and venous systems interact via a third porous medium, representing "small scale" vasculature.This compartment is given a simplified representation in terms of only a throughput conductivity K c (x) and a porosity ϕ c (x), leading to fluxes and transit times given by Δτ c = ϕ c (x)/q(x), i.e., pore volume per throughput volume rate.While subscript c alludes to capillaries, the compartment will typically also capture other parts of the microvascular network like arterioles and venules.Although the capillary compartment per construction excludes any lateral transport, by regulating the throughput conductivity it may induce such transport in the arterial or venous compartments.
For model identification purposes, we consider eight parameter fields, three scalar porosity fields ϕ a , ϕ v and ϕ c , two anisotropic diagonal permeability tensor fields K a and K v , and one scalar permeability field K c .(For a 3D case there will be ten parameter fields.)The spatial domain is divided into a given number of regions of interest, N ROI , which gives a total number of 8 � N ROI parameters to estimate (10 � N ROI in 3D).Observation of time-series for the contrast agent concentration on each voxel will guide the identification process.As outlined above, the local perfusion estimate for each voxel is taken to be the capillary flow rate q.
The spatial domain will be embedded in a regular Cartesian mesh, where selected cells will be made inactive to map out more complex geometrical structures.Boundary conditions will be a combination of pressure conditions and no-flow conditions, and the flow pattern will typically be driven by the difference between an arterial inlet pressure and a venous outlet pressure.Using a five (seven in 3D) point finite volume stencil, we solve for pressures p a and p v , velocities u a and u v , and transfer flux q.
To solve for the transport of contrast agent, we utilize that the formulation chosen for the fluid flow obeys a discrete maximum principle.Thus any two voxels connected by a streamline have a strict upstream-downstream relation and the contrast agent advection can be computed sequentially voxel by voxel starting from the most upstream cell.We combine mass conservation r � u = q and time-of-flight τ, defined by u � rτ = ϕ into r � (τ u) = ϕ + τq, and consider a bounded control region of volume V. We partition the mass rate, F, through the region boundaries according to the flow entering (up) and leaving (dn): Thus one attribute an average arrival time, τ up for fluid entering, and an average departure time, τ dn for fluid leaving.Assuming the source term q to be uniform over the volume, the associated average time τ q will either carry an upstream value (q > 0) into a venular compartment, or leave (q < 0) an arterial compartment with a value "somewhere"between τ up and τ dn .For the computation in the current paper, we have lumped the transit times into a single quantity for each compartment, and write Dt v ¼ V� v =F dn v for q > 0 (venular) and Dt a ¼ V� a =F up a for q < 0 (arterial).This essentially means that we consider the flow to enter capillaries at the downstream end of an arterial compartment, and depart capillaries at the upstream end of the venular compartment.The contrast agent concentration is sampled at discrete times of uniform increment Δt, and the advection of the discrete contrast vector is carried out in three steps related to the arterial, capillary, and venous compartment respectively.Starting from the most upstream arterial voxel, the upstream concentration c up a ðtÞ is shifted by transit time Δτ a to obtain the downstream concentration c dn a ðtÞ ¼ c up a ðt À Dt a Þ.We also compute volume-averaged concentration c vol a ðtÞ, serving the dual purpose of arterial voxel concentration and input concentration for the capillary compartment.For multidimensional domains, the upstream concentration c up a ðtÞ will generally be a flux-weighted combination over several upstream neighbor connections, confer the lumping process leading to Eq 2. Similarly, by shifting and averaging using the capillary transit time Δτ c , we obtain the capillary downstream concentration

A hybrid-scale flow model used to generate synthetic data
Here we consider a hybrid-scale model formulation, integrating one-dimensional network flow with a two-or three-dimensional continuum model.Known ground truth parameter values are assumed for this model when generating synthetic data.We presented the full details of this model in Hodneland et.al., [15], and here we will briefly recap some main points.Compared with the pure Darcy formulation presented in the preceding section, the arterial and venous compartments are now both augmented with flow networks representing visible arterial and venous vessel structure respectively.The nodes of each network are classified as interior, terminal, and root nodes.A single node can have multiple roles.We associate a pressure value to each node location, and neighboring nodes are connected by edges representing vessel segments.For each edge, we then have a one-dimensional flux that is related to node pressures via the Poiseuille law.We impose a mass balance condition at each node, requiring associated fluxes to satisfy net zero accumulation.
Each root node represents a cut-off for our region of interest, and we assign pressure boundary conditions.Each terminal node represents a coupling to the associated continuum model, and we considered the flux to be proportional to the pressure difference between the node and the continuum at each node location.In order to provide for a smooth transition, a distribution function with finite support is centered at each location.We chose a distribution function with finite support under the assumption that individual terminals cannot substantially provide large tissue slabs with blood but are rather local in nature.However, this assumption suggests that the behavior of the tracer may become abrupt and discontinuous at certain points where the support ends, which might not accurately represent the real-world scenario and must be considered an approximation.In particular, for small support radii, i.e. small values of � in [15,Eq 15], the use of a finite support function may lead to sharp fronts of the tracer.The same approach is pursued in [16].An alternative solution is to use a smooth function like a Gaussian with infinite support and gradual decay, which means it can provide a more realistic representation of transitions and gradients, but still be unphysiological in terms of maximum range.
We combine the above with standard Darcy formulations for the arterial and venous continuums respectively, and introduce a coupling term between the two in terms of a flux proportional to local pressure difference.We identify this flux as the perfusion supplying the local tissue.We solve the total coupled formulation to obtain pressure and flux distribution across the combined model.Based on the fluxes obtained, we model the advection of the contrast agent.At each root node of the arterial network, we specify contrast concentration in the entering bloodstream in terms of an arterial inflow function.At each terminal node, we again apply the Gaussian kernel to obtain a smooth distribution of the entering contrast agent.Compared with the previous section where we used an explicit approach for agent transport, we here perform a fully implicit procedure based on a backward Euler discretization of the advection equation.

Bayesian inversion
The Bayesian method we use to assimilate tracer concentration data is introduced in [11].This methodology has proven to be robust and efficient when assimilating big amounts of data to large-scale sub-surface reservoir models, see e.g., [17,18].The technique is based on a regularized Levenberg-Marquardt method and iteratively searches for the minimum of an average cost problem.The approach utilizes an ensemble of model parameters, where each member of the ensemble is denoted m j , j = 1, . .., N, and N is the ensemble size.The unknown parameters are tissue properties such as transmissibility and porosity.Further, we denote the blood circulation model M; R N m !R N d and the simulated concentration data are denoted c j ¼ Mðm j Þ; j ¼ 1; . . .; N. The vector c j contains concentration data at all points in time where measurements are available and for all voxels in the discretized domain.In each voxel the concentration is given by c vol total (see the 'Methods' section).The cost function of interest is given by: where c o is the tracer concentration data and m b j are given background (or best guess) parameters (defined below).The second term in the cost function is the regularization part, and γ is a (user-defined) weight parameter.The data errors are assumed to be Gaussian-distributed with zero mean and covariance C c .The matrix Cm is the sample covariance of the background parameters.The cost function can be derived from the theory for Bayesian inference, which states that the posterior distribution for a random variable, given knowledge of data, is proportional to the prior probability distribution times the likelihood distribution, i.e., p(m|d) / p(m) p(d|m).Given assumptions about Gaussianity for the distributions, the first term in the cost function comes from the likelihood part and the second (regularization) term is a consequence of the prior part of Bayes' formula, and prevents overfitting of (noisy) data.For more details, we refer to [19,20].
The above problem is solved iteratively and, without going into details (see [11] for a full method description), we state the update formula for ensemble member j at iteration i + 1: where m i j � m b j is adaptively updated.Initially, these parameters are randomly drawn from a user-defined distribution (later referred to as the initial ensemble).The columns in the matrices S i m and S i c are given by m i j À m i and c i j À Mðm i Þ, respectively, where m i is the sample mean of the parameter vectors at iteration i.The perturbation terms ϵ j are samples from the data error distribution.Similar to the background parameters, the weighting term γ i is updated for each iteration.This parameter is increased if the average data mismatch increases, otherwise the parameter is reduced.Generally, large values for γ i result in smaller steps that are suitable for highly nonlinear problems, and small values for γ i implies larger steps that are suitable to problems that are closer to linear.In our application the values are selected based on experience with this specific problem, but also based on experience with the estimation of parameters in Darcy models in general.In addition, a stopping criteria is selected such that the data mismatch of the final solution is of the same magnitude as the noise level.
Localization.As mentioned in the introduction, inaccurate correlations between parameters and measurements make it necessary to restrict the updates of parameters.The inaccurate correlations are due to sampling errors that comes from the limited ensemble size we use to estimate correlation matrices.For example, concentration data in a particular voxel is likely to have a low correlation with tissue properties in voxels that are far from that data position.Since sampling errors in the correlation matrices that are estimated as part of Eq 4 are inevitable, localization will limit unrealistic updates far from the data positions.We pursue a correlation-based localization technique [13], which is best described by formulating the update step (Eq 4) as where where p < N is the number of singular values after truncation.The matrix U p 2 R N d �p contains the left-singular vectors, W p is a diagonal matrix composed of the singular values, and V p 2 R N�p contains the right-singular vectors.We keep 99.9% of the sum of descending singular values.Substituting Eq 8 in Eq 5 gives T Dc i j is the projection of the difference between simulated and real measurements, onto a subspace consisting of the dominant directions contained in U p .The matrix Ãi is the Kalman gain matrix after data projection.Introducing a tapering matrix Λ i , we define the update step where � denotes the matrix Schur product.The tapering matrix is either 1 or 0, and is (at a given iteration i) computed based on information about correlations between the projected data innovations (Dc i j ) and the model parameters (m i j ).The details for computing the tapering matrix are not included here, and we refer to [13] for a full description.Examples of localization areas are given below.

Workflow
We summarize the methodological description with a schematic workflow, Fig 1 .First, an initial ensemble must be specified, based on prior knowledge about the tissue or organ that is investigated.The ensemble is usually randomly sampled from statistical probability distributions.Next, the iterative ensemble smoother is run, where execution of the forward circulation model and parameter updates are done alternately.The update part requires measurements as input, which is generated synthetically in our application.Finally, after the iterative smoother has converged and posterior parameters are obtained, the perfusion is returned from a final execution of the circulation model.

Results
In the current work, we illustrate the methodology on a two-dimensional problem.The example is based on an image of a frog tongue found in the book by Cohnheim [21], see Fig 2 .The frog tongue is stretched and pinned to the surface at six points so that the thickness is sufficiently small for visual identification of major arteries and veins.The dimension of the domain we have set to be 33 mm in the x-direction, and 40 mm in the y-direction.This data has the advantage of being truly 2D but still biologically relevant.Selecting a 2D slice from a 3D geometry will not be sufficient to account for out-of-plane fluxes.The measurements are synthetically generated using the hybrid-scale model, and the unknown parameters are estimated using the porous media model.The code and data are available at [22].
With frequent sampling of data in time and space, the amount of data would be excessive.Therefore, only a fraction of the data is used.Since time-correlations for data errors are hard to accurately specify (although for a different scientific research area, a useful discussion regarding time-correlated data errors can be found in [23]), an algorithm for finding a reasonably reduced set of data is designed.The algorithm is designed to use a fraction of the points in time (i.e., it uses whole images for the selected time points).In the case study described below, we decided to use one-tenth of the available time points.In our case, the hybrid-scale model is sampling data every second for a time interval of 150 seconds.To select the time points we form one vector x i for each time point which we consider as a stochastic variable.Then we form the covariance matrix C = Cov(x i , x j ) where 1 � i, j � n t and n t = 150 denote the number of time points where data are sampled.Let I be a set of indexes from the set {1, . .., n t } and let C I,I denote the principal submatrix of C containing the rows and columns from the set I. Motivated by D-optimal design [24] we would like to select n p = bn t /10c points such that det C I,I is maximized over all sets I having n p distinct entries from the set {1, . .., n t }. (For clarification: b�c denotes the floor operation, 10 is our choice, and might obviously be changed.) When n t is increasing, the maximization problem described above would be very time-consuming to solve by a brute-force search.Since solving this optimization is only done to improve the speed and efficiency of the Bayesian inversion, we have designed a greedy approach to find a reasonable solution to the problem of finding the set I which maximizes det C I,I over all sets I of size n p .The greedy algorithm is designed by first selecting k 1 such that C k 1 ;k 1 is maximized.This will be the solution to the maximization problem if the size of the set I is 1.From this, we start to form sets I 0 of increasing size until we reach the size n p by adding one index in each step.In the second step, we find k 2 with k 1 fixed from the first step such that det C I 0 ,I 0 where I 0 has size 2. For an arbitrary step we are given indexes k 1 , . .., k j−1 and find k j maximizing det C I 0 ,I 0 where I 0 = {k 1 , . .., k j−1 } [ k where k is any of the indexes 1 � k � n t different from the previous selected indexes k 1 , . .., k j−1 .The algorithm terminates when the size of I 0 is n p .
The experimental setup ("truth") for the hybrid-scale model is given in Table 1.The model is run on a grid with dimension 515 × 634, and data are contrast agent concentration values in every grid cell.As a base case, the simulated measurements are upscaled using a stencil of 4 × 4 grid cells.In order to limit the amount of data we down-sample the measurements as described above.This gives a total of 191220 values.The resulting spatial-temporal resolution is within what can be achieved with both MRI and CT scans [25].The bloodstream entering the domain (AIF) is represented by a Gamma-variate function, see Fig 3.
In this example study we select to estimate log-transformed transmissibilities (transmissibilities are defined by T = K/μ) and porosities in the arterial, venous, and capillary compartments.Although scalar transmissibilities are used for the "true" model we do not want to assume that this information is known in general, and in the arterial and venous compartments we estimate transmissibility in both the x-and y-direction.This gives a total of eight unknown parameter fields.Initial values are generated as Gaussian random fields with given mean and Gaussian variograms.We assume constant mean values in the frog tissue, and have selected values with moderate differences from the "true" properties listed in Table 1, except for the capillary conductivity (K c ) which is substantially larger.The rationale for this choice Table 1.Experimental setup for the hybrid-scale model.See [15,26] for more information about the model.stems from the fact that this quantity is directly related to the perfusion, and we want to demonstrate the robustness of the methodology in a case with limited prior knowledge of the capillary (or throughput) conductivity.One might argue that there exist accurate prior knowledge if the methodology is applied to a human organ, but the ability to adapt to a large variation of values for capillary conductivity is still of high importance in order to capture the properties of damaged tissue caused by e.g., stroke or cancer.In addition, there are variations between individuals caused by cardiovascular health, age, lifestyle factors, underlying medical conditions, medications, etc., that makes it difficult to provide accurate prior properties.The standard deviations for the distributions reflect the uncertainty of the prior information.In our study, we select 0.1 for porosity values and 1 for log-transformed transmissibilities.The capillary porosity (ϕ c ) determines the transit time through the capillaries.This quantity is not present in the hybrid model, and for simplicity, we have selected a low prior value and low standard deviation.For real data, the prior distributions should be tailored to each specific application, and based on the knowledge of the organ or tissue that is studied.The selection of prior distributions is not crucial for finding better estimates for the parameters, but some discrepancies for the updated models are reported [20] The variogram ranges are drawn from Gaussian distributions with mean values equal to 74 and 60 grid blocks in the x-and y-direction, respectively.The variance (in both directions) is equal to one grid block.The relatively long variogram ranges are based on the assumption that the tissue has an approximately homogenous structure.

Parameter
In the arteries and veins, we use transmissibility (in both x-and y-direction) given by T = D 2 /32μ, where D is the vessel diameter.This relation is found by comparing the Darcy equation (with porosity equal to 1) with the Hagen-Poiseuille equation for laminar incompressible flow.
To ensure that all parameters are physically reasonable we also impose upper and lower bounds for the values.The bounds are specified rather loosely in order to capture spatial variations of tissue properties.We have approximated the tissue with a porous media model, which means that vessels of different diameters are represented using different values for permeability and porosity in the model.
The statistical parameters are listed in Table 2.For completeness, we have included the permeabilities (in parenthesis) in addition to the log-transformed transmissibilities in the table.The porous media flow model is discretized using 128 cells in the x-direction, and 158 cells in the y-direction (base case).There are 7507 inactive cells that fall outside the boundary of the frog tongue.Each parameter field then consists of 12717 values, and the total number of parameters is N m = 101736.We have also listed the true parameters (see also Table 1) used to generate the measurements with the hybrid-scale model.The viscosity, boundary pressures and AIF are identical for the hybrid-scale model and the porous media model.
On Fig 4 we show three realizations from the the prior distribution for y-transmissibility in the arterial compartment, and the porosity in the venous compartment.Note that the blood vessels are not visible on the porosity fields, because the mean porosity values are equal (0.1) both inside and outside the vessels.This value is adopted because it is used in the hybrid-scale model when generating the measurements, and the motivation is to avoid porosity values larger than 1 in the parts of the tissue where arteries and veins are crossing (on top of each other).Examples of the computed localization domains based on the initial ensemble, l 0 kl , for selected observations, are shown on Fig 5 .The areas change for each projected observation (l), and at every iteration (i).The initial ensemble is used as starting values for the algorithm given by Eq 4, and we use 200 ensemble members (model realizations).The maximum number of iterations is set to 10, and in addition, the algorithm will stop if the relative change in average data mismatch is less than 10%.We assume that the measurement error standard deviation for the data is 10% of the maximum observed concentration and we assume that the noise is uncorrelated in time and space.This means that C c is a diagonal matrix.The initial value for the weight parameter γ 0 is 1, and the reduction and increment factors are 0.9 and 2, respectively.These factors are used to reduce or increase γ i at each iteration.In this example, five iterations were always performed before the method terminated.

Estimated perfusion
The estimated perfusion follows from Eq 1.The true perfusion (see [15], Equation 7) is shown on the left picture on Fig 6 .The perfusion range is between 5 mL/min/100mL and 21 mL/min/ 100mL.The spatial relative error in estimated perfusion is shown in the middle and right pictures.The relative error is computed as (P T − P E )/P T , where P T is the true perfusion and P E is the estimated perfusion.The prior and posterior errors are calculated by simulating the perfusion using the mean of the prior and posterior ensembles.Note that the perfusion is zero in the vessels, and the vessels are colored for enhanced visibility.There is a clear improvement for the posterior estimate, and the relative error is especially reduced in the lower part of the domain, where the prior values are as low as −14%.In particular, large errors occur when the true perfusion is low.The relative error using the posterior estimate is minimum −2.6%.From the results, we see that both the prior and posterior overestimate the perfusion in most of the domain.In order to further quantify the accuracy of the estimated perfusion we compare the results with two widely used techniques within radiology: the maximum slope method and deconvolution based on singular value decomposition.We will not present these methods here, but refer to [27] for readers interested in details.In Fig 7 (left) we show the average perfusion (P) computed for the entire frog tongue domain.The true value is shown in blue and is 9.5 mL/min/100mL.The maximum slope method and deconvolution method slightly overestimate the perfusion value (the error is 1.1 mL/min/100mL), whereas the prior model returns a value as high as 71 mL/min/100mL.The posterior model obtained after the assimilation of concentration data is 8.9 mL/min/100mL, slightly below the true value.On the right plot, we show corresponding results for four regions, obtained by dividing the domain into equally where P q T is the true perfusion in region q and P q E is the estimated perfusion in the region q, using one of the estimation (E) techniques: maximum slope, deconvolution, prior model, or posterior model.In the above calculations, we utilize the information about the position of arteries and veins by setting the perfusion value in a region to zero if the region contains more than 50% of vessels.The maximum number of regions are all voxels in the base case domain (128 × 158), excluding the voxels that contain arteries or veins.The conclusion from this figure is that the posterior model is better than all the other techniques we have considered, for all partition levels.
The computational time (measured as wall clock seconds) for the different grids used in the above computations are shown in Table 3.The forward circulation model simulations of the ensemble realizations are not parallelized in this study and it is therefore a potential to speed up the simulations significantly.For this study, the calculations take less than 25 minutes to   finish for the most refined grid.Preliminary studies on a 3D domain with grid 45 × 32 × 35 (number of regions equal to 50400) shows that one iteration (with an ensemble size equal to 200) takes approximately 20 minutes.However, the forward simulations dominate the time to finish one iteration and the benefit of parallelizing the simulations increase rapidly when the gird discretization is refined and the number of regions increase.

Estimated concentration
where the measurements (c o ) and simulated observations (c i j ) are concentration data.

Discussion
We have demonstrated the use of ensemble-based estimation techniques for assimilating processed contrast agent concentration data from dynamic MRI acquisitions into models for blood perfusion in organs.The methodology is applied to a domain with known vascular structure, and the realism of the methodology is enhanced by selecting different models for synthetic data generation and data assimilation.The data-generating hybrid-scale model includes a detailed representation of both the arterial and venous vessel structure, and a distribution function is used for the flux from the vessel network to the continuum.The porous media model used for data assimilation is simpler and needs to represent vessel structures as high-permeable channels.We show that improved perfusion estimates are obtained, compared to traditional methods (maximum slope and deconvolution methods).The major contribution of our work is the ability to provide better estimates of the perfusion in any region of interest embedded in the organ of investigation.The importance of accurate estimates of perfusion is, as mentioned in the introduction, wide-ranging.The work presented here represents a first step towards a full-scale tool for clinical usage.The next step is to evaluate the methodology on synthetic three-dimensional contrast data, and eventually real MRI-based measurements.When using real measurements there are challenges related to how the concentration data are scaled and pre-processed, that must be addressed.In a real setting, there is also uncertainty in the boundary conditions, i.e., the boundary pressures and AIF.A natural extension of the methodology is to augment the set of unknown parameters with these quantities and update the boundary conditions as part of the assimilation workflow.It has recently been suggested that adding in-silico trials would be helpful for the development of better treatment of acute ischemic stroke [28][29][30] and new modeling approaches to handling different issues are being developed [31].Obviously, patient-specific models that align with observations would increase the value of in-silico models even further.
c dn c ðtÞ ¼ c vol a ðt À Dt c Þ and the capillary volume average c vol c ðtÞ.For the venular compartment, we again start from the most upstream voxel and for each voxel, we compute a lumped upstream concentration vector c up v ðtÞ as a flux weighted sum over upstream neighbor connections and the out-flux c dn c ðtÞ from the capillary compartment.Shifting by transit time Δτ v yields the downstream concentration vector c dn v ðtÞ ¼ c up v ðt À Dt v Þ. Finally the venous volume average concentration c vol v ðtÞ is computed, and combined into a total (including extravascular tissue) volume concentration for each voxel c vol total ðtÞ ¼ � a c vol a ðtÞ þ � c c vol c ðtÞ þ � v c vol v ðtÞ.
and I N d is the N d -dimensional identity matrix.Correlation-based localization involves the computation of a correlation matrix between observations and model parameters.In order to avoid huge memory requirements, we use a truncated singular value decomposition (TSVD) to project the data onto a subspace consisting of the dominant singular vectors.Applying TSVD to Si c gives Si

Fig 2 .
Fig 2. Vascular network for a stretched frog tongue.Arteries are red, and veins are blue.The gray color represents the tongue tissue.The inlets and outlets are at the lower edge.https://doi.org/10.1371/journal.pcbi.1011127.g002

Fig 8 .
Fig 8. Mean absolute error for estimated perfusion for increasing number of regions.On the plot, 'MS' are the maximum slope results, 'DC' are the deconvolution results, 'PR' are the results using the prior model, and 'PO' are results using the posterior model after the data assimilation.The x-axis is logarithmic, and the range is from 1 to 10483.https://doi.org/10.1371/journal.pcbi.1011127.g008

Fig 9 .
Fig 9. Spatial distribution of contrast concentration (mmol/L) at four points in time.Blue color indicates low concentration and yellow indicates high concentration.The same scale (colormap) is used for all plots.The columns show from left to right the data, the prior mean, the posterior mean, and the posterior standard deviation.https://doi.org/10.1371/journal.pcbi.1011127.g009

Fig 9
Fig 9 shows the spatial distribution of contrast concentration at four points in time, simulated using the full resolution with 128 × 158 grid cells.At time equal to 33 seconds a clear concentration pattern is seen in the data.The prior estimate illustrates the fact that the blood vessels in the porous media model are not sealed, and there is diffusion through the vessel walls.This phenomenon is reduced in the posterior model, but the circular concentration areas seen in the data are not recovered.The rightmost figure shows the standard deviation for the posterior ensemble.As explained in the introduction, it is preferable to avoid underestimation of the uncertainty after iterating.The standard deviations are higher than zero, especially early in the time period, thereby indicating that collapse of the ensemble is avoided.Visually, we obtain improvements for the posterior concentration at all points in time, compared to the prior concentration values.See also S1 Video.The estimated concentration is further visualized at nine positions in the frog tongue, see Fig 10.These positions are evenly distributed in the domain, and among the points, the middle right is in a high permeability zone (blood vessels).Fig 11 shows concentration curves at the nine positions.A clear reduction in the uncertainty is seen for the posterior curves (shown in green), compared to the prior curves (shown in blue).The

Fig 11 .
Fig 11.Temporal distribution of contrast concentrations at nine positions, given by the (coarse) grid indices in the titles.Red color indicates the measurements, blue color represents the initial ensemble, and green color represent the updated posterior ensemble.The x-axis is the time in seconds, and the y-axis is the concentration values in mmol/L.https://doi.org/10.1371/journal.pcbi.1011127.g011

Table 2 . Statistical parameters used to generate the initial ensemble.
Lower and upper bounds are denoted lb and ub, respectively.Note that the mean transmissibility values are only valid in the frog tissue (outside the identified main vessels).Inside the vessels, the transmissibilities are given by the Hagen-Poiseuille equation.The values in parentheses are permeabilities corresponding to the log-transformed transmissibilities.