Viscosity quantification using multi-contrast magnetic particle imaging

The measurement signal in MPI are voltages U(t) induced into one or more receive coils at times $t\in [0,T)$ where T is the period length of the applied measurement sequence. For simplicity we consider a single receive coil in this work. In a pre-processing step, the periodic time signal U is commonly expanded into a Fourier series with coefficients ${\hat{u}}_{k}={\int }_{0}^{T}U(t){{\rm{e}}}^{-2\pi {\rm{i}}{k}/T}{\rm{d}}t$, which allows one to filter out the feed through of the excitation signal and noisy frequency components prior to improve reconstruction [23].

The standard assumption made in MPI is that the spatial distribution of the non-interacting particles $c(\overrightarrow{r})$ and the frequency components of the measured signal can be described by an integral equation linear in the particle distribution

$$\begin{eqnarray}&&{\hat{u}}_{k}={\int }_{{\rm{\Omega }}}{\hat{s}}_{k}(\overrightarrow{r})c(\overrightarrow{r}){{\rm{d}}}^{3}r,\end{eqnarray} \tag{ 1 }$$

where ${\rm{\Omega }}\subset {{\mathbb{R}}}^{3}$ is the spatial domain and ${\hat{s}}_{k}(\overrightarrow{r})$ is the system function [1, 21]. We note that this linear assumption is experimentally proven to be valid, if the particle concentration is sufficiently low [24].

Beside linearity, one key assumption made in (1) is that the magnetic response of the particles can be described by a single system function. However, this assumption is easily violated, e.g. by temporal or spatial changes of the particle temperature, viscosity or binding state, or by having multiple tracers differing in their particle core size or particle hydrodynamic size distributions inside the field of view.

As shown in [3] the standard model can be too restrictive in cases where multiple particle types are used for imaging. We therefore propose to incorporate these cases by introducing an additional parameter space $V\subseteq {{\mathbb{R}}}^{l}$, $l\in {\mathbb{N}}$ to describe all experimentally relevant nanoparticle states such that each parameter $\overrightarrow{\tau }\in V$ is in a one-to-one correspondence to a nanoparticle state. Experimentally relevant parameters are those possibly changing the MPI signal such that in a scenario where e.g.the temperature is constant throughout an experiment, one may very well omit the temperature from the parameter space. I.e. V is a subset of a l-dimensional vector space, where each dimension encodes a degree of freedom such as the temperature, the viscosity of the particles environment, different particle types or microscopic details such as the diameter of the particle core. To shorten notation, we will denote particles in a state described by the parameter $\overrightarrow{\tau }\in V$ as $\overrightarrow{\tau }$-particles.

This additional parameter space then extends the spatial domain to ${{\mathbb{R}}}^{3}\times V$ on which a generalized particle distribution ${\mathfrak{c}}(\overrightarrow{r},\overrightarrow{\tau })$ is defined, which describes the density of $\overrightarrow{\tau }$-particles at each location $\overrightarrow{r}\in {{\mathbb{R}}}^{3}$. The generalized particle distribution is linked to the particle distribution of the single-contrast signal equation (1) by $c(\overrightarrow{r})={\int }_{V}{\mathfrak{c}}(\overrightarrow{r},\overrightarrow{\tau }){{\rm{d}}}^{l}\tau$. Likewise, a corresponding generalized system function ${\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},\overrightarrow{\tau })$ can be defined, which describes the magnetic response of $\overrightarrow{\tau }$-particles. The generalized signal equation is then given by

$$\begin{eqnarray}&&{\hat{u}}_{k}={\int }_{{\rm{\Omega }}}{\int }_{V}{\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},\overrightarrow{\tau }){\mathfrak{c}}(\overrightarrow{r},\overrightarrow{\tau }){{\rm{d}}}^{l}\tau {{\rm{d}}}^{3}r.\end{eqnarray} \tag{ 2 }$$

We note that (2) is formulated in a very general manner. Often the generalized particle distribution ${\mathfrak{c}}$ will have additional properties, which change the mathematical structure of equation (2). Three cases where a property lead to a simplification of the generalized imaging equation are discussed next.

2.2.1. Separable parameter dependency

Consider microscopic parameters such as the core size distribution of the particle core and the hydrodynamic particle diameter. Assuming that this parameter does not change with the spatial location of the particles the generalized particle distribution

$$\begin{eqnarray}&&{\mathfrak{c}}(\overrightarrow{r},\overrightarrow{\tau })=c(\overrightarrow{r})\varrho (\overrightarrow{\tau })\end{eqnarray} \tag{ 3 }$$

factorizes into the spatial particle distribution c and a state density ϱ describing the microscopic details, with ${\int }_{V}\varrho (\overrightarrow{\tau }){{\rm{d}}}^{l}\tau =1$. In this case a system function ${\hat{s}}_{k}(\overrightarrow{r}):= {\int }_{V}{\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},\overrightarrow{\tau })\varrho (\overrightarrow{\tau }){{\rm{d}}}^{l}\tau$ incorporating all the microscopic details can be defined. Inserting ${\mathfrak{c}}$ into equation (2) and using the definition of ${\hat{s}}_{k}(\overrightarrow{r})$ one obtains the MPI standard model formulated in equation (1) as a special case of our generalized model. Since most particles used in MPI have a polydisperse core size and hydrodynamic size distribution, (3) is a common assumption.

2.2.2. Discrete set of parameters

Second, the parameter space is assumed to be discrete. I.e. there is only a finite number of experimentally relevant states within the parameter space, e.g. immobilized particles marking a medical instrument and a mobile tracer within the blood pool [25] or different particle types [3]. Let ${\overrightarrow{\tau }}_{i}\in V$, i = 1, 2, ..., n be the parameters describing the given particle types and ${c}_{{\overrightarrow{\tau }}_{i}}$ be the corresponding distribution of ${\overrightarrow{\tau }}_{i}$-particles. In this case one can link the generalized particle distribution ${\mathfrak{c}}$ with the individual distributions ${c}_{{\overrightarrow{\tau }}_{i}}$ by

$$\begin{eqnarray}&&{\mathfrak{c}}(\overrightarrow{r},\overrightarrow{\tau })=\displaystyle \sum _{i}{c}_{{\overrightarrow{\tau }}_{i}}(\overrightarrow{r})\delta (\overrightarrow{\tau }-{\overrightarrow{\tau }}_{i}).\end{eqnarray} \tag{ 4 }$$

Within this discrete setting one can define a corresponding system function ${\hat{s}}_{k,{\overrightarrow{\tau }}_{i}}(\overrightarrow{r})=\hat{{\mathfrak{s}}}(\overrightarrow{r},{\overrightarrow{\tau }}_{i})$ for each ${\overrightarrow{\tau }}_{i}$ such that the signal equation (2) simplifies to

$$\begin{eqnarray}&&{\hat{u}}_{k}=\displaystyle \sum _{i}{u}_{k,{\overrightarrow{\tau }}_{i}}=\displaystyle \sum _{i}{\int }_{{\rm{\Omega }}}{\hat{s}}_{k,{\overrightarrow{\tau }}_{i}}(\overrightarrow{r}){c}_{{\overrightarrow{\tau }}_{i}}(\overrightarrow{r}){{\rm{d}}}^{3}r.\end{eqnarray} \tag{ 5 }$$

In this model that was first considered in [3] for multi-contrast imaging, the measurement signal ${\hat{u}}_{k}$ is the superposition of a finite number of individual signals ${u}_{k,{\overrightarrow{\tau }}_{i}}$ induced by the particles of type ${\overrightarrow{\tau }}_{i}$ as described by the MPI standard model.

2.2.3. Field-like parameters

Last, one can consider field-like properties like e.g. temperature and viscosity, which take a specific value at each spatial position. In this case the state parameter is coupled to the spatial position via a vector-valued function $\varphi \,:{\rm{\Omega }}\to V$. I.e. all particles at $\overrightarrow{r}\in {\rm{\Omega }}$ are in the state $\overrightarrow{\tau }=\varphi (\overrightarrow{r})$. Therefore, ${\mathfrak{c}}(\overrightarrow{r},\overrightarrow{\tau })$ is zero for all other states and equal to the local particle density $c(\overrightarrow{r})$ for $\overrightarrow{\tau }=\varphi (\overrightarrow{r})$. We can express this in a distributional sense by

$$\begin{eqnarray}&&\hat{{\mathfrak{c}}}(\overrightarrow{r},\overrightarrow{\tau })=c(\overrightarrow{r})\delta (\overrightarrow{\tau }-\varphi (\overrightarrow{r})).\end{eqnarray} \tag{ 6 }$$

Provided the map φ is fix, the signal generated by a particle distribution c is then given by

$$\begin{eqnarray}&&{\hat{u}}_{k}={\int }_{{\rm{\Omega }}}{\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},\varphi (\overrightarrow{r}))c(\overrightarrow{r}){{\rm{d}}}^{3}r,\end{eqnarray} \tag{ 7 }$$

where ${\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},\varphi (\overrightarrow{r}))$ could be seen as a conventional system matrix.

Image reconstruction in MPI aims to invert the forward problem given by the MPI signal equation. It involves the solution of an inverse problem by appropriate methods. In general one can differentiate between algebraic methods that treat the signal equation in a general way [21, 26, 27] and analytic methods that exploit its structure and derive direct reconstruction formulas [5, 28]. A good overview about these different approaches can be found in [29]. In the current work we formulate the reconstruction formulas in an algebraic sense although it would also be possible to apply analytical methods, like the x-space reconstruction.

The goal of a single-contrast image reconstruction is to obtain $c(\overrightarrow{r})$ from measurements ${\hat{u}}_{k}$ by inverting equation (1). Similarly, the goal of a multi-contrast image reconstruction is to determine ${\mathfrak{c}}(\overrightarrow{r},\overrightarrow{\tau })$ from the measurements ${\hat{u}}_{k}$ by inverting equation (2). However, since ${\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},\overrightarrow{\tau })$ is usually nonlinear in $\overrightarrow{\tau }$ we expect the related inverse problem to be even more difficult to tackle than the single-contrast image reconstruction problem. Even if a solution to the inverse problem exists, iterative nonlinear solvers are not guaranteed to find it since the solver might get stuck at local minima. Prior to discussing the most general case, the single-contrast reconstruction will be recapped and strategies for handling the multi-contrast signal equation for discrete parameter sets (5) and field-like parameters (7) will be discussed.

2.3.1. Single-contrast image reconstruction

As initially proposed in [1], MPI reconstruction was performed in a single-contrast fashion by solving the inverse problem (1) for a given measurement and system function. In practice, the signal equation is first transformed into a discrete linear system of equations, which is then solved by linear least squares minimization [1, 21, 22]. From the point of view of our generalized signal equation (2) single-contrast image reconstruction uses a fixed system function ${\hat{s}}_{k,{\overrightarrow{\tau }}_{m}}(\overrightarrow{r})={\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},{\overrightarrow{\tau }}_{m})$ and optimizes the distribution ${c}_{{\overrightarrow{\tau }}_{m}}(\overrightarrow{r})={\mathfrak{c}}(\overrightarrow{r},{\overrightarrow{\tau }}_{m})$ of ${\overrightarrow{\tau }}_{m}$-particles,

$$\begin{eqnarray}&&{c}_{{\overrightarrow{\tau }}_{m}}=\mathop{\mathrm{argmin}}\limits_{\tilde{c}}\displaystyle \sum _{k=1}^{M}| {u}_{k}-{\int }_{{\rm{\Omega }}}{\hat{s}}_{k,{\overrightarrow{\tau }}_{m}}(\overrightarrow{r})\tilde{c}(\overrightarrow{r}){{\rm{d}}}^{3}r{| }^{2}.\end{eqnarray} \tag{ 8 }$$

Prior to solving this linear system of equation it is reduced by filtering out noisy frequency components. The reduced least squares problem can then be efficiently solved by regularized and weighted iterative methods as discussed in [27].

However, despite the common usage of (8) it is not ensured that the least squares minimizer ${c}_{{\overrightarrow{\tau }}_{m}}(\overrightarrow{r})$ equals the actual distribution of particles $c(\overrightarrow{r})$. This is due to a possible mismatch between ${\overrightarrow{\tau }}_{m}$ that was used during the calibration of ${\hat{s}}_{k{\overrightarrow{\tau }}_{m}}$ and the parameter $\overrightarrow{\tau }$ that is valid during the object measurement. In a sense, ${c}_{{\overrightarrow{\tau }}_{m}}(\overrightarrow{r})$ provides the best-approximation of the measurement signal when considering a distribution of ${\overrightarrow{\tau }}_{m}$-particles only. It has already experimentally been proven that a mismatch between the system function $\overrightarrow{\tau }$ and the object $\overrightarrow{\tau }$ leads to serious artifacts and non-quantitative results [30].

2.3.2. Multi-channel image reconstruction for discrete parameter sets

Consider a setting with a discrete parameter set as described by signal equation (5), i.e. the measurement signal ${\hat{u}}_{k}$ can be written as a superposition of the individual signals of a finite number of MPI tracers. In this case the single-contrast reconstruction approach can be straight forwardly extended by taking a finite number of calibrated system functions ${\{{\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},{\overrightarrow{\tau }}_{m})\}}_{m=1}^{l}$ into account [3]. The optimization problem (8) can then be reformulated to

$$\begin{eqnarray}&&{c}_{{\overrightarrow{\tau }}_{1}},\,...,{c}_{{\overrightarrow{\tau }}_{l}}=\mathop{\mathrm{argmin}}\limits_{{\tilde{c}}_{{\overrightarrow{\tau }}_{1}},\ldots ,{\tilde{c}}_{{\overrightarrow{\tau }}_{l}}}\displaystyle \sum _{k=1}^{M}| {u}_{k}-\displaystyle \sum _{m=1}^{l}{\int }_{{\rm{\Omega }}}{\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},{\overrightarrow{\tau }}_{m}){\tilde{c}}_{{\overrightarrow{\tau }}_{m}}(\overrightarrow{r}){{\rm{d}}}^{3}r{| }^{2}.\end{eqnarray} \tag{ 9 }$$

In case the calibration measurements ${\hat{s}}_{k{\overrightarrow{\tau }}_{i}}={\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},{\overrightarrow{\tau }}_{i})$ $\forall i$ match the particles during the measurements, the reconstructed distributions ${c}_{{\overrightarrow{\tau }}_{m}}$ equal the respective spatial distribution of ${\overrightarrow{\tau }}_{m}$-particles. In case the MPI tracers used during the experiment are not matched with the system functions, one cannot directly link the generalized distribution ${c}_{{\overrightarrow{\tau }}_{m}}$ to a spatial particle distribution. One major challenge of the multi-contrast reconstruction method as formulated in (9) is that an increase of the number of system functions (i.e. discretization steps l) will have a negative influence on the ill-posedness of the inverse problem.

2.3.3. Multi-channel image reconstruction for field-like continuous parameter spaces

The signal of measurement scenarios with field-like continuous parameter spaces is described by signal equation (7). In this scenario image reconstruction aims to reconstruct both, the distribution of the particles and the vector-valued function φ describing the spatial dependence of the parameter $\overrightarrow{\tau }\in V$

$$\begin{eqnarray}&&c,\varphi =\mathop{\mathrm{argmin}}\limits_{\tilde{c},\tilde{\varphi }}\displaystyle \sum _{k=1}^{M}| {u}_{k}-{\int }_{{\rm{\Omega }}}{\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},\tilde{\varphi }(\overrightarrow{r})){\tilde{c}}_{l}(\overrightarrow{r}){{\rm{d}}}^{3}r{| }^{2}.\end{eqnarray} \tag{ 10 }$$

Since $\overrightarrow{\tau }$ directly describes the properties of the tracer it makes no sense assign $\overrightarrow{\tau }$ some specific value at a location where there is no tracer. Thus one would impose further restrictions upon φ to ensure the uniqueness of a solution.

The main issue with the optimization (10) is that ${\hat{{\mathfrak{s}}}}_{k}$ needs to be known. Since measuring ${\hat{{\mathfrak{s}}}}_{k}$ for a single parameter $\overrightarrow{\tau }$ is already quite tedious, a measurement based approach as for the single-contrast scenario seems to be unpractical. Alternatively, one could simulate ${\hat{{\mathfrak{s}}}}_{k}$ using a model based approach. However, current models are either too simple to reproduce the rich magnetization dynamics of the particles or numerically too costly to use.

2.3.4. General multi-channel image reconstruction

Likewise, the aim of the general multi-channel image reconstruction is to obtain the generalized distribution ${\mathfrak{c}}(\overrightarrow{r},\overrightarrow{\tau })$ from a given measurement u_k and generalized system function ${\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},\overrightarrow{\tau })$

$$\begin{eqnarray}&&{\mathfrak{c}}=\mathop{\mathrm{argmin}}\limits_{\tilde{{\mathfrak{c}}}}\displaystyle \sum _{k=1}^{M}| {u}_{k}-{\int }_{{\rm{\Omega }}}{\int }_{V}{\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},\overrightarrow{\tau })\tilde{{\mathfrak{c}}}(\overrightarrow{r},\overrightarrow{\tau }){{\rm{d}}}^{l}\tau {{\rm{d}}}^{3}r{| }^{2}.\end{eqnarray} \tag{ 11 }$$

Here too the main problem is the current unavailability of a suitable generalized system function.

As done in [20], one can circumvent this problem by a discretization of the parameter space

$$\begin{eqnarray}&&{\int }_{{\rm{\Omega }}}{\int }_{V}{\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},\overrightarrow{\tau }){\mathfrak{c}}(\overrightarrow{r},\overrightarrow{\tau }){{\rm{d}}}^{l}\tau {{\rm{d}}}^{3}r\approx \displaystyle \sum _{m=1}^{l}{\int }_{{\rm{\Omega }}}{\hat{{\mathfrak{s}}}}_{k}(\overrightarrow{r},{\overrightarrow{\tau }}_{m}){c}_{{\overrightarrow{\tau }}_{m}}(\overrightarrow{r}){{\rm{d}}}^{3}r,\end{eqnarray} \tag{ 12 }$$

where ${c}_{{\overrightarrow{\tau }}_{m}}(\overrightarrow{r})={w}_{m}{\mathfrak{c}}(\overrightarrow{r},{\overrightarrow{\tau }}_{m})$ with appropriate discretization weights w_m. This discretization reduces the general image reconstruction problem in equation (11) to the multi-channel image reconstruction for discrete parameter sets (5), yielding l reconstructed distributions ${c}_{{\overrightarrow{\tau }}_{m}}$. For instance in [20] two measured system functions for particles of low temperature ${\overrightarrow{\tau }}_{\mathrm{low}}$ and high temperature ${\overrightarrow{\tau }}_{\mathrm{high}}$ were used to obtain temperature and particle distribution from the corresponding reconstruction channels ${c}_{{\overrightarrow{\tau }}_{\mathrm{low}}}$ and ${c}_{{\overrightarrow{\tau }}_{\mathrm{high}}}$.

Within the experimental part of this work we will proceed quite similarly. We will obtain two system functions at low and high viscosity and perform the multi-channel reconstruction in equation (11). The mapping between the signal distribution of the reconstructed channels and the viscosity will be obtained by a calibration procedure. A detailed description of this procedure can be found in section 4.

A series of viscosity samples was prepared by mixing ferucarbotran (Resovist, FUJIFILM RI Pharma Co., Ltd, Japan) with an iron concentration of 500 mmol_Fe l⁻¹, glycerol and distilled water. The first seven samples M1 to M7 have log spaced viscosities 1.0 mPa s, 1.93 mPa s, 3.73 mPa s, 7.20 mPa s, 13.9 mPa s, 26.8 mPa s, 51.8 mPa s, whereas the last three U1, U2, and U3 have viscosities of 2.8 mPa s, 20.4 mPa s, 100 mPa s at a room temperature of 23 °C. A total amount of V = 50 μl was prepared per mixture, each with an iron concentration of c_Fe= 50 mmol_Fe l⁻¹.

As there is no adequate theoretical description on the viscosity of water-glycerol mixtures we use an empirical formula describing the viscosity of a water-glycerol mixture η [31]. It is given by

$$\begin{eqnarray}&&\eta ={\eta }_{w}^{\alpha }{\eta }_{g}^{1-\alpha },\end{eqnarray} \tag{ 13 }$$

where η_w and η_g are the respective viscosities of water and glycerol and α is a weighting factor depending on the temperature and the glycerol mass fraction ${C}_{{m}_{g}}$. This formula was numerically inverted to obtain ${C}_{{m}_{g}}$ for each of the desired mixture viscosities. At a lab temperature of 22 °C the densities of water and glycerol are ϱ_w = 997.6 kg m⁻³ and ϱ_g = 1263 kg m⁻³. The glycerol volume fraction of the mixture is given by

$$\begin{eqnarray}&&{C}_{{V}_{g}}=\displaystyle \frac{{\varrho }_{w}{C}_{{m}_{g}}}{{\varrho }_{w}{C}_{{m}_{g}}+{\varrho }_{g}(1-{C}_{{m}_{g}})},\end{eqnarray} \tag{ 14 }$$

where ${C}_{{m}_{g}}$ is the glycerol mass fraction and the corresponding mass and volume fraction of water are given by ${C}_{{m}_{w}}=1-{C}_{{m}_{g}}$ and ${C}_{{V}_{w}}=1-{C}_{{V}_{g}}$. With these the samples were prepared in three steps.

• 1.  
  Ferucarbotran was diluted with distilled water to lower the iron concentration down to $\tfrac{{c}_{\mathrm{Fe}}}{{C}_{{V}_{w}}}$ to ensure that the final mixture had the desired iron concentration of c_Fe, where the assumption is made that the diluted ferucarbotran has the same viscosity as distilled water. Next, a volume of $V(1-{C}_{{V}_{g}})$ was pipetted from the diluted ferucarbotran.
• 2.  
  A mass of ${{VC}}_{{V}_{g}}{\varrho }_{g}$ of glycerol was weighted using a high precision balance instead of having to pipette the extremely viscous liquid.
• 3.  
  Diluted ferucarbotran and glycerol were mixed and filled into a glass capillary with an inner diameter of 1.3 mm. The total height of the mixture within the capillary was about 9.4 mm.

An overview over all samples with their viscosity, glycerol mass fraction, glycerol volume fraction and iron concentration of the ferucarbotran part can be found in table 1.

Table 1.  The table shows a complete list of all 50 μl samples prepared with their desired viscosity η. Additionally, the glycerol mass and volume fraction ${C}_{{m}_{g}}$ and ${C}_{{V}_{g}}$, the iron concentration c_Fe and volume of the water part of the mixture and the mass of the glycerol part required for the preparation of the samples are listed.

Sample name	M1	M2	M3	M4	M5	M6	M7	U1	U2	U3
η(mPa s)	1.0	1.93	3.73	7.20	13.9	26.8	51.8	2.8	20.4	100
${C}_{{m}_{g}}$	0.018 2	0.248 4	0.417 6	0.546 3	0.647 6	0.729	0.796 2	0.349 9	0.697 2	0.852 3
${C}_{{V}_{g}}$	0.014 4	0.207	0.362	0.488	0.592	0.680	0.755	0.298	0.645	0.820
cFe (mmolFe l−1)	50.7	63.1	78.3	97.6	122.6	156.3	204.4	71.3	140.9	278.0
Vw (μl)	49.3	39.6	31.9	25.6	20.4	16.0	12.2	35.1	17.7	9.0
mg (mg)	0.910	13.071	22.831	30.781	37.382	42.932	47.685	18.836	40.735	51.777

In this paragraph the relative uncertainty of the viscosity will be calculated using the law of propagation of errors. This uncertainty is mainly influenced by uncertainties during the sample preparation, temperature changes during the experiment, and uncertainties due to the empirical nature of equation (13). Whenever the propagated error of a source depends on the viscosity the maximum relative error is taken as a bound for all samples.

The preparation of the desired glycerol mass fraction is mainly influenced by the accuracy of the pipette, the accuracy of the balance and temperature deviations during sample preparation. Pipetting could be done with a relative error of no more than 0.7% for the smallest volumes of 11.6 μl. The balance used had a relative error of 0.01%. The pipetting error directly propagates to a relative 0.7% error of the glycerol mass fraction ${C}_{{m}_{g}}$ and the weighting step accumulates an relative error of 0.01%. Temperature uncertainties slightly affect the pipetting step due to changes in the water density. However, even uncertainties of 4 °C would only lead to a relative error of about 0.1% in ${C}_{{m}_{g}}$. Since all three sources of error are independent from each other the total relative error of the desired glycerol mass fraction is 0.7%. The propagated relative viscosity error stemming from the glycerol mass fraction depends on the value of ${C}_{{m}_{g}}$ and can be bound by 4%.

The temperature deviations inside the scanner bore were recorded using a fiber optical thermometer (Fotemp, Weidmann Technologies Deutschland GmbH, Germany). During measurements temperature ranged from 21.5 °C to 22.5 °C. These changes cause a uncertainty in the weighting factor α and the viscosities of glycerol and water. Propagation of the 0.5 °C uncertainty yields a relative error of 6%. According to [31] the average relative prediction error of the viscosity formula (13) is below 2.3%. Accumulation of these independent errors yields a total relative uncertainty in the viscosity of 7.5%.

MPI measurements were carried out using a preclinical MPI scanner (Philips preclinical MPI package with a Bruker preclinical MPI system). The system was operated with a gradient field of 1.0 T μ₀⁻¹ m⁻¹, 1.0 T μ₀⁻¹ m⁻¹, $2.0\,{\rm{T}}\,{\mu }_{0}^{-1}\,{{\rm{m}}}^{-1}$ in x-, y-, and z-direction and a drive-field amplitude of 12 mT μ₀⁻¹ in x- and z-direction, respectively, where μ₀ is the vacuum permeability. In our system the x-direction coincides with the direction of the bore, the z-axis points upwards and the y-direction is the remaining direction in the horizontal plane orthogonal to the bore.

Particle excitation is done by moving the field free point on a 2D Lissajous trajectory using excitation frequencies of 2.5/102 MHz and 2.5/99 MHz. This results in a cycle duration of 652.8 μs and a drive-field field of view of 24 × 12 mm in the xz-plane. During each measurement two coils in x- and z-direction pick up the induced voltage signal. These two signals are then filtered and digitized with a bandwidth of 2.5 MHz. The final measurement vector is a concatenation of the Fourier transformed time signals. For each measurement a sample was attached to the arm of a three axis robot (Isel Automation GmbH) with the glass capillary orthogonal to the excitation plane to ensure precise placement within the field of view.

In total two system matrix and ten regular measurements were performed. For general information on MPI measurements we refer to [32]. The two calibration measurements were obtained on a 26 × 12 mm system matrix field of view (SM-FoV) with a voxel size of 1 × 1 mm by averaging 5000 MPI measurements at each grid position for the samples M1 and M7, respectively. For each of the two samples the measurements were collected and post processed into a system matrix. Henceforth we will refer to the system matrix of sample M1 by S₁ and the one of M7 by S₇. For the regular measurements each of the samples was placed in the center of the SM-FoV, obtaining 5000 MPI measurements for each one.

From the set of all recorded frequencies only those above 80 kHz and with a signal to noise ratio above 5 are selected. The remaining frequencies are removed from the measurement vectors and system matrices. To assess the influence of noisy data the 5000 measurements were block averaged with a block length M_avg of 100, 50, 20, 10, 5, 2, 1 resulting in a total of 50, 100, 250, 500, 1000, 2500, 5000 reconstructed multi-channel images per sample.

Using the two system matrices S₁ and S₇ a multi-channel image reconstruction proposed in [3] was performed. I.e. for each measurement vector $\hat{{\bf{u}}}$ the Tikhonov regularized optimization problem

$$\begin{eqnarray}{{\bf{c}}}_{1},{{\bf{c}}}_{7}=\mathop{\mathrm{argmin}}\limits_{{\tilde{{\bf{c}}}}_{1},{\tilde{{\bf{c}}}}_{7}\in {{\mathbb{R}}}_{+}^{N}}{\parallel \left(\begin{array}{cc}{{\bf{S}}}_{1} & {{\bf{S}}}_{7}\end{array}\right)\left(\begin{array}{c}{\tilde{{\bf{c}}}}_{1}\\ {\tilde{{\bf{c}}}}_{7}\end{array}\right)-\hat{{\bf{u}}}\parallel }_{2}^{2}+\tilde{\lambda }{\parallel \left(\begin{array}{c}{\tilde{{\bf{c}}}}_{1}\\ {\tilde{{\bf{c}}}}_{7}\end{array}\right)\parallel }_{2}^{2}\end{eqnarray} \tag{ 15 }$$

is solved. Here, N = 26 × 12 × 1 is the total number of voxels within each channel, $\tilde{\lambda }$ is the regularization parameter and c₁ and c₇ are the channels of the multi-channel reconstruction. Usually the regularization parameter is not reported directly, but the parameter $\tilde{\lambda }=\lambda N/\parallel \left(\begin{array}{cc}{{\bf{S}}}_{1} & {{\bf{S}}}_{7}\end{array}\right){\parallel }_{{\rm{F}}}$, where $\parallel \cdot {\parallel }_{{\rm{F}}}$ denotes the Frobenius-norm. Equation (15) is solved iteratively with 10000 iterations of the Kaczmarz solver and $\tilde{\lambda }={10}^{-5}$.

Each reconstructed image contains two channels c₁ and c₇ as shown in figure 1, with the sample signal located in the center. For each image the signal from each channel is collected by summation over a 6 × 4 mm region of interest in the center. Note that the width of the ROI is larger than its height to account for the lower gradient strength in x-direction. We refer to the summed signals from channel c₁ as σ₁ and to channel c₇ as σ₇, respectively. Depending on the averaging block size M_avg one ends up with 50, 100, 250, 500, 1000, 2500, 5000 values σ₁ and σ₇ for each sample, as shown in figure 2.

In the last processing step the channel weight

$$\begin{eqnarray}&&\lambda ({\sigma }_{1},{\sigma }_{7})=\displaystyle \frac{{\sigma }_{1}}{{\sigma }_{1}+{\sigma }_{7}},\end{eqnarray} \tag{ 16 }$$

is calculated. It assigns each tuple (σ₁, σ₇) a number in the interval $[0,1]$ which is invariant under rescaling $({\sigma }_{1},{\sigma }_{7})\mapsto (\alpha {\sigma }_{1},\alpha {\sigma }_{7})$, $\alpha \in {{\mathbb{R}}}_{+}$. The idea is to quantify the relative signal weight of the channels. With λ = 1 only c₁ contains signal from a sample, whereas with λ = 0 only c₇ contains signal. This definition is purely heuristically and motivated by the observation that the relative signal weight contained within the channels shifts from c₁ to c₇ as the viscosity of the sample increases as shown in figure 1. The complete mapping from MPI measurement data to $\lambda \in [0,1]$ is summarized in figure 3.

Zoom In Zoom Out Reset image size

Viscosity estimation is based on a calibration curve. The curve is obtained by using every fourth data point from the data sets of M1 to M7 with M_avg = 100 as training data for a piece wise linear curve. All remaining data, including ${M}_{\mathrm{avg}}\ne 100$ is used for testing of the viscosity estimation. The corresponding values of λ(σ₁, σ₇) of the training data are averaged for each sample and paired with the samples viscosity. As shown in figure 4 the values λ decrease monotonously with increasing viscosity.

Based on the data points we define a monotonous function λ(η) using linear interpolation in between and linear extrapolation beyond the data points as shown in figure 4. This in turn allows one to define the inverse mapping η(λ) as a simple root finding problem

$$\begin{eqnarray}&&\eta (\lambda ({\sigma }_{1},{\sigma }_{7})):= \tilde{\eta }\ \mathrm{such}\ \mathrm{that}\ \lambda (\tilde{\eta })=\lambda ({\sigma }_{1},{\sigma }_{7}).\end{eqnarray} \tag{ 17 }$$

This mapping is used to estimate the viscosities of the remaining data sets as shown exemplary in figure 5. For each sample and block averaging number the mean viscosity η_estim, the standard deviation σ and systematic deviation ${\rm{\Delta }}\eta =| {\eta }_{\mathrm{sample}}-{\eta }_{\mathrm{estim}}|$ from the sample viscosity η_sample were calculated as listed in table 1.

All multi-channel reconstructions yield a two-channel image with the signal from the sample located in the center as shown in figure 1. We observed a strong correlation between block averaging length M_avg and image noise, with low noise for M_avg = 100 and increasing noise for smaller M_avg. Moreover, the multi-channel reconstruction maps signal from low viscosity samples M1, M2, and U1 predominantly to the channel M1 and signal from high viscosity samples M6, M7, and U3 predominantly to the channel M7. For samples with medium viscosity M3, M4, M5, and U2 one observes quite significant changes in the signal distribution between both channels as shown in figure 1. Though only one reconstruction is shown, this trend can be observed for all reconstructed data sets. In most cases this allows one to visually discriminate low, medium, and high viscosity samples. Due to improved signal to noise ratio discrimination works better for large block averaging length. Due to the fact that the signal shifts nonlinearly between the channels, we were not able to estimate the viscosity solely based on the visual perception of the reconstructed two-channel images.

Zoom In Zoom Out Reset image size

Further processing of each image yields signal weights σ₁ and σ₇. Due to the aforementioned noise in the reconstructed images one observes a more or less wide spread distribution of σ₁ and σ₇. The distribution is wider for low M_avg and narrower for large M_avg as shown in figure 2 for each sample, M_avg = 1 and M_avg = 100. Independent of the sample one observes a reduction of the distribution width by a factor of about seven as the block averaging number increases from M_avg = 1 to M_avg = 100. The only exception is sample M1 with M_avg = 50 and M_avg = 100, where the width of the distribution is exactly zero. As λ(σ₁, σ₇) is calculated directly from σ₁ and σ₇ we find that the observations above also apply for the relative signal weight.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next we study the mean value of σ₁, σ₇, and λ(σ₁, σ₇). These mean values mainly depend on the sample and show only slight variations with the block averaging length. For the mean value of σ₁, σ₇, we observe a non-monotonous change with respect to the sample viscosity as shown in figure 2, where the samples are sorted with respect to their viscosity. A quantitative analysis reveals that on average σ₁ is larger for sample M2 than for sample M1. From there on σ₁ decreases as the viscosity increases. Likewise, the mean signal of σ₇ increases monotonously from sample M1 to M6 and decreases from there on, i.e. for the last two samples M7 and U3. Contrary, the mean relative signal weight λ(σ₁, σ₇) decreases with increasing viscosity for block averaging lengths of 5, 10, 20, 50, 100, which is the basis of our viscosity estimation.

Using the method proposed in this work, we are able to estimate the viscosity of all samples. Due to the noise in the signal the estimated viscosities are distributed around a mean value for all block averaging length and samples. The distributions are wider for smaller M_avg and narrower for larger M_avg as shown in figure 5 for block averaging length of 1, 10, 100. The mean estimated viscosity, its standard deviation and its systematic deviation from the sample viscosity are listed in table 2.

Zoom In Zoom Out Reset image size

A closer analysis reveals that the viscosities of the samples M1 to M7, which were also used in the calibration process, were best approximated. On average estimated and sample viscosity differed up to 3%, 2%, 3%, 5%, 12%, 40%, 54% for block averaging length of 100, 50, 20, 10, 5, 2, 1, respectively. Down to M_avg = 10 the sample viscosity lies within an interval of one standard deviation around the mean estimated viscosity. Assuming independence of systematic and statistical errors a total methodological error of the estimation of up to 6%, 8%, 10%, 15%, 21%, 50%, 64% for block averaging length of 100, 50, 20, 10, 5, 2, 1 can be calculated.

As for the samples U1, U2, and U3, which were not used for calibration we need to distinguish between the samples U1 and U2 with viscosities within the range of the calibration samples and U3 with a viscosity outside the calibration rage. For the viscosity estimation of U1 and U2 one observes a higher average difference between estimated and sample viscosity for block averaging length of 100, 50, 20, 10 and lower ones for block averaging length of 5, 2, 1. In total the methodological errors are up to 6%, 8%, 10%, 13%, 18%, 25%, 28% block averaging length of 100, 50, 20, 10, 5, 2, 1. The viscosity of sample U3 is significantly underestimated with its prepared value being 100 mPa s and its estimated value ranging between 62 mPa s for M_avg = 100 to 30 mPa s for M_avg = 1. Performing a separate analysis, we found that an inclusion of U3 into the set of calibration samples reduces this problem to the point, where the sample U3 is as good approximated as U1 and U2.