{'fraction_non_alphanumeric': 0.03863170389554332, 'fraction_numerical': 0.013650588108341426, 'mean_word_length': 4.713625154130703, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Evaluation of three inverse problem models to quantify skin microcirculation using diffusion-weighted MRI

Skin microcirculation plays an important role in diseases such as chronic venous insufficiency and diabetes. Magnetic resonance imaging (MRI) can provide quantitative information with a better penetration depth than other noninvasive methods, such as laser Doppler flowmetry or optical coherence tomography. Moreover, successful MRI skin studies have recently been reported. In this article, we investigate three potential inverse models to quantify skin microcirculation using diffusion-weighted MRI (DWI), also known as q-space MRI. The model parameters are estimated based on nonlinear least-squares (NLS). For each of the three models, an optimal DWI sampling scheme is proposed based on D-optimality in order to minimize the size of the confidence region of the NLS estimates and thus the effect of the experimental noise inherent to DWI. The resulting covariance matrices of the NLS estimates are predicted by asymptotic normality and compared to the ones computed by Monte-Carlo simulations. Our numerical results demonstrate the effectiveness of the proposed models and corresponding DWI sampling schemes as compared to conventional approaches.


Motivation
The skin is known as the largest organ of the human body, both by weight and surface area, and plays a critical role in protecting the body against pathogens found in the environment. It is comprised of three primary layers: the epidermis (from 0.05 mm on the eyelids to 1.5 mm on the palms and soles), the dermis (from 0.3 mm on the eyelids and 3.0 mm on the back), and the hypodermis (subcutaneous adipose layer), from the outermost to the innermost. Blood vessels in the dermis provide nourishment and waste removal for its own cells, as well as the cells in the deepest layers of the epidermis (Stratum basale) via diffusion. These blood vessels are composed of capillaries, arterioles, venules, and arteriovenous anastomosis (shunting vessels). The hypodermis is used mainly for fat storage and contains larger blood vessels and nerves.
Perfusion through capillaries pertains to nutrition, while perfusion through the other types of blood vessels refers to temperature regulation and feeding and draining of the capillary network [1]. Skin perfusion can be impaired by diseases such as chronic venous insufficiency [2] and diabetes [3], which can lead to decubitus ulcer formations and necrosis. In the wound healing process, vascular supply to the wound is essential and relies on neovascularization or angiogenesis [4], and means of quantitatively characterizing perfusion is crucial for prognosis.

Experimental methods for skin microcirculation
Laser Doppler flowmetry (LDF) is one of the most widely used noninvasive methods to monitor microvascular blood flow, even though LDF is limited to a sampling depth of about 1 mm due to light scattering and LDF measurements are non-absolute [1]. Recent efforts to remedy the latter require the a priori knowledge of the geometry and optical properties of the microvasculature [1], and yet do not contain any information about flow directionality [5].
Optical coherence tomography (OCT) is a noninvasive optical technique that can provide morphological information with an isotropic spatial resolution of about 15 µm in a 10×10 mm 2 field-of-view (FOV) and a penetration depth limited to 0.5 to 1.5 mm due to light scattering. Thus, OCT can only reach down to the superficial layer of the dermis [6,7]. Confocal microscopy can achieve spatial resolutions of about 1 µm and resolve the different layers of the epidermis at the cellular level [7]. When combined with confocal Raman spectroscopy, additional biochemical information targeting specific skin structures becomes available [8]. Fluorescence fiber-optic confocal microscopy requires the injection of a fluorophore, and can produce in vivo images of epidermis cell layers down to the subpapillary dermis with subcellular resolution [9].
Recent advances in hardware now allow the use of magnetic resonance imaging (MRI) and spectroscopy (MRS) for in vivo skin studies on state-of-the-art clinical whole-body scanners. Specific surface coil can provide higher signal-to-noise ratios (SNR) than volume coils, allowing an elementary volume (voxel) on the order of 20 × 100 × 800 µm 3 for imaging, and a volume of interest of 140 × 1,000 × 1,000 µm 3 for spectroscopy [7,10]. Unlike LDF, MRI can be quantitative and the FOV can be controlled. Moreover, diffusion-weighted MRI (DWI) can make use of the diffusion of water molecules in biological tissues to probe microstructures which cannot be seen using standard imaging. DWI protocols involve the sampling of q-space (the Fourier reciprocal space of spin displacements) via the choice of magnetic gradient strengths and directions [11,12,13], while the imaging is still done in k-space (the Fourier reciprocal space of the spin locations). Several review articles and book chapters describe the uses of DWI (e.g., [14,15]). In particular, DWI methods can provide quantitative information about the displacement spectrum of water molecules and average velocities, such as the intra voxel incoherent motion method (IVIM [16]). We are proposing developing DWI protocols specifically for the study of blood flow in the skin microvasculature. Since the relaxation times of skin tissue are relatively low compared to other tissues [10,17], stimulated-echo pulse sequences may be preferred to spin-echo sequences to avoid an excessive drop of signal due to T * 2 relaxation [18,19]. The minimal voxel size allowable depends on the achievable SNR, which in turn varies according to the hardware available. Preliminary experiments will be required to determine this voxel size, which is expected to be on the order of 200×800×800 µm 3 , i.e., somewhere between the minimum resolution achievable for standard imaging and spectroscopy.
This manuscript aims at evaluating three potential models that are based on volume-averaged transport equations for blood flow in the skin vasculature and would be used in conjunction with DWI protocols to extract quantitative information about skin microcirculation.

Transport equation for volume-averaged skin microcirculation
The transport of a solute of concentration C(x, t) in a medium is governed by the advectiondiffusion equation. By assuming incompressible flow, the standard advection-diffusion equation can be volume-averaged and linearized [20] for each MRI voxel, resulting in is the spatially averaged velocity field, and D ∈ R 3×3 is the symmetric dispersion tensor that reflects on both the advection and diffusion processes, as well as on the geometry, and may be time-dependent [21]. V and D are both defined on a voxel-by-voxel basis.
Let us now consider specific models for skin microcirculation. The imaging voxel can be taken such that z is the direction perpendicular to the skin surface (and surface coil) and is assumed to be one of the principle axis of the dispersion tensor. As a result, dispersion along the z direction is decoupled from dispersion in the (x, y)-plane, such that Moreover, as we attempt to characterize steady or quasi-steady conditions, the net z velocity component should vanish (i.e., V z ≡ 0), provided the FOV contains the outermost capillary layer. Otherwise, blood would either be leaving or accumulating in the blood vessels in the z-direction. Therefore, the governing equation for C(x, t) can be written as It can be shown (e.g. [22]) that the response to an impulse function, C(x, y, z, 0) = C 0 δ(x, y, z, 0) after a time t is given by For an impulse initially located at a location x 0 , x in equation (4) needs to be replaced by the vector x − x 0 , which can be thought of as a displacement vector.

DWI protocol for skin microcirculation
Pulsed field gradient imaging results in the acquisition of a normalized echo attenuation, E(q, ∆) ∈ C, which corresponds to the Fourier transform of the probability density function of the spin displacements [11]. Computed as the Fourier transform of C(x, t) with C 0 = 1, E(q, ∆) is a function of the spatial frequency of the spin displacements, q := γgδ/2π ∈ R 3 , and ∆, where δ and ∆ are the duration of and time between the pulsed field gradients g ∈ R 3 , respectively [11]. By decomposing q into its component in the (x, y)-plane, with In practice, E(q, ∆) is obtained by dividing the measured echo with diffusion-encoding gradients turned on (q = 0) by the echo obtained with zero diffusion gradients (q = 0). The characterization of blood microcirculation in the skin is obtained by fitting the experimental normalized echo attenuation to equation (5), which constitutes an inverse problem. Blood flow is characterized via six parameters: the four coefficients of the the dispersion tensor (D xx , D yy , D xy , D zz ), and the two average velocity components parallel to the skin surface (V x , V y ). The physical interpretation of the dispersion tensor measured via q-space DWI is slightly different from the dispersion terminology introduced earlier. In engineering processes, only dispersion within flow ducts is considered, while the MRI signal originates from both the blood vessels and the surrounding tissue, since the latter also contains protons. Equation (5) may therefore need to be completed by a contribution due to the surrounding tissue. The uncertainty lies in the T 1 -and T 2 -weighting factor, which may render this contribution negligible and depends on the imaging pulse sequence. Accounting for this additional contribution, equation (5) becomes where f blood and f surr are the T 1 -and T 2 -weighted volume fractions ("visible" volume fractions as in [19]) of blood and surrounding tissue within the voxel, respectively, with f blood + f surr = 1. E blood (q, ∆) is given by (5), and E surr (q, ∆) is assumed to take the form of an anisotropic unrestricted diffusion compartment, and D surr is a diffusion tensor of the form given in (2) with four unknown coefficients based on symmetry arguments (z is a principle direction). Three models are then proposed for the study of skin microcirculation using MRI protocols with increasing levels of complexity: The echo attenuation is given by Eq. (5)-(6), which contains six (6) unknown parameters to determine. Either the surrounding does not contribute to the MRI signal or can be lumped into the dispersion tensor description. Model 2. The echo attenuation is given by Eq. (7)-(8) with D surr reduced to D surr I, corresponding to isotropic diffusion inside the cutaneous tissue surrounding the microvasculature. Model 2 contains eight (8) unknown parameters to determine. Model 3. The echo attenuation is given by Eq. (7)-(8) with D surr taking the form given in (2), which contains eleven (11) unknown parameters to determine. The surrounding contributes to the MRI signal as an anisotropic unrestricted diffusion compartment with the same symmetry as the microvasculature.
The visible volume fractions, f m , are equal to the actual volume fractions, only if the relaxation times of the different diffusion compartments are the same, which is what was explicitly assumed in [16]. Otherwise, relaxation times can be estimated based on additional or prior measurements (e.g., [23,17]) in order to compute the "true" volume fractions after the extraction of the visible volume fractions, f m , using DWI data [19].

Optimal experiment design
The objective of this manuscript is then to numerically find the optimal q-space sampling scheme for the inverse problems stated in section 2.2, and demonstrate the feasibility and efficiency of such models subject to realistic noise levels. Noise in MRI images is typically modeled as independent Gaussian noise with zero mean in the two quadrature detection channels [24,25]. The sensitivity to noise will determine the required SNR in experimental DW images, and thus, ultimately, the voxel size. This numerical analysis of our models will therefore allow the proper design of MRI protocols for future experiments, and avoid costly trial-and-error repetitions of MRI experiments. Consider the true model given as E(q m , ∆, θ ⋆ ), with the true parameter vector θ ⋆ ∈ R p .  Here, E is the expectation operator and * denotes the complex conjugate operator. The origin of this complex-valued white noise is from the asymptotic normal distribution of the Fourier transform of white noise (see more details in [26,27]). It may also come from the quantization and electronic noise of the data acquisition system.
Given the measurements (9), we consider to find the nonlinear least squares (NLS) estimatê θ M that minimizes the residual sum of squares: where | · | denotes the magnitude operator on a complex number. In relating Q(M ) to the covariance matrix of the NLS estimates, it can be shown that the Fisher Information Matrix (FIM) of the model (9) is given by [28]: where Re{·} operates element-wise and E ′ (q m , ∆, θ ⋆ ) := ∂E(q m , ∆, θ)/∂θ| θ=θ ⋆ . By the Cramer-Rao theorem [29], the covariance matrix of any unbiased estimatorθ is lower bounded by the Cramer-Rao Lower Bound (CRLB), or the inverse of the Fisher Information Matrix FIM(θ ⋆ , Q): Moreover, with some regularity conditions [28,30],θ M is asymptotically normal as M → ∞, θ M → d N (θ ⋆ , W ), with mean θ ⋆ and covariance matrix where Σ := lim M →∞ Σ 2M and Σ is a positive definite matrix. In practice, the exact θ ⋆ is not available to evaluate Σ 2M in (11). We only have θ 0 ≃ θ ⋆ as an a priori knowledge. Also, on clinical scanners, because of the limited gradient strength capabilities, the values for the duration (δ) and separation (∆) of the pulsed gradients are very constrained. The maximum value for q is typically achieved by using the maximum allowed gradient strength, which then fixes δ. The value for ∆ is then typically chosen in order to minimize the duration of the echo time in order to minimize T 2 relaxation. Therefore, the problem is to find a correct configuration of M measurement points, Q := {q m } M m=1 such that the size of the confidence region of the parameter estimates is minimized. This is achieved by maximizing the determinant of the FIM (D-optimality [29]): This optimization problem is solved numerically by using the optimization toolbox in Matlab for a given θ ⋆ . In practice, the a priori value of θ ⋆ may be used instead of the unknown true value for θ ⋆ [28].

Results
Our three models for skin microcirculation are evaluated by analyzing the covariance matrices from the NLS estimation using Monte-Carlo (MC) simulations with K realizations and a noise level σ = 5%. The parameter estimation is done using the Levenberg-Marquardt algorithm [31] with an initial guess equal to a vector with a 10% deviation from the true parameter θ ⋆ . The number K of MC realizations is chosen such that the predicted standard deviations converge towards their predicted values. As shown in figure 1(a), the normalized error has converged after K ≈ 5000 realizations. Indeed, for K = 5000 the error in the trace of W is below 0.7% for model 1 (below 1.8% for model 2, and below 0.4% for model 3).   For each model, two sets of DWI protocols identified by sets of locations on one or more spheres are compared: our optimized DWI sampling scheme (Q opt ) with a given number of points on each sampling sphere, and a conventional sampling scheme based on Papadakis' scheme 12 [32] distributed over each sphere (Q P ), as was done experimentally in [19] for instance. We used one sphere for model 1, four spheres for model 2 and six spheres for model 3, resulting in different numbers of optimized sampling locations (M = 12, 48, and 72, respectively). These DWI sampling schemes are considered to be realistic since they would result in reasonable acquisition times for a patient as only a single slice at a low resolution is acquired. For Q opt , the gradient magnitudes and distribution over the predetermined number of spheres were optimized according to (14). For Q P , unlike [19], the gradients magnitudes were optimized for each sphere. The comparison bears on the mean values (Table 1) and the standard deviations (Tables 2 and 3) for the NLS parameter estimation based on covariance matrices and analytical predictions using the asymptotic normality of the NLS estimates (13).  Our first observation is that the NLS estimates of the model 1 are consistent. For model 2 and 3, the correct mean values are obtained within a tolerable bias error given the finite number of sampling locations (Table 1). In terms of parameter estimation uncertainty, Table 2 gives the comparison of the three models according to the standard deviation for Q opt . The addition of a second compartment in models 2 and 3 penalizes the standard deviation of all the parameters. Nevertheless, the MC estimation stays very close to the prediction concerning the dispersion tensor, with less than 4.5% difference for model 2. Results using Q P are shown in Table 3. Standard deviations of NLS estimates resulting from the optimized schemes are smaller than those from the conventional schemes, with improvements up to 40%. For further comparison of our proposed models, we are interested in knowing how many sampling locations are needed for models 2 and 3 in order to have similar standard deviations as for model 1. Based on our predictions, since the convergence of MC has been proven, we can plot the trace of the prediction matrix as a function of the number of measurements. On figure 1(b), we can see that only very large number of sampling locations can allow models 2 and 3 to have the same standard deviation as model 1, which is not feasible in practice. However, a 90-point scheme for model 3 yields the same uncertainty as a 48-point scheme for model 2.

Conclusions
Three models based on a linearized volume-averaged advection-diffusion equation are proposed to quantify skin microcirculation using DWI. For each model, the DWI protocols are optimized based on D-optimality in order to counteract the inherently low SNR caused by the short relaxation times of skin tissues and reduce the uncertainty in the NLS parameter estimates. Our models and our optimized vs. conventional DWI protocols are validated via MC simulations. This manuscript is thus meant as a tool to design quantitative DWI protocols for skin microcirculation studies (e.g., to investigate chronic venous insufficiency or diabetes) as a complement to the ongoing hardware advances.