Multimodality characterization of microstructure by the combination of diffusion NMR and time-domain diffuse optical data

Combining datasets with a model of the underlying physics prior to mapping of tissue provides a novel approach improving the estimation of parameters. We demonstrate this approach by merging near infrared diffuse optical signal data with diffusion NMR data to inform a model describing the microstructure of a sample. The study is conducted on a homogeneous emulsion of oil in a dispersion medium of water and proteins. The use of a protein based background, rich in collagen, introduces a similarity to real tissues compared to other models such as intralipids. The sample is investigated with the two modalities separately. Then, the two datasets are used to inform a combined model, and to estimate the size of the microstructural elements and the volume fraction. The combined model fits the microstructural properties by minimizing the difference between experimental and modelled data. The experimental results are validated with confocal laser scanning microscopy. The final results demonstrate that the combined model provides improved estimates of microstructural parameters compared to either individual model alone.


Introduction
Non-destructive investigation of material microstructure is a general challenge for imaging science. Microstructure imaging exploits the sensitivity of the signal obtained from an imaging device to the microscopic architecture of materials, such as emulsion droplet sizes (van Staveren et al 1991, Denkova et al 2004, material inclusion in minerals (Kaiser et al 1997) or indication of biological properties (Schmitt and Kumar 1998), to estimate features of the microstructure within a sample. In biomedical imaging these ideas lead to the possibility of non-invasive histology. Histology provides the gold standard for diagnosis of a wide range of diseases (Schonberg et al 2006), but is invasive and limited to small target areas. Microstructure imaging potentially offers similar information, non-invasively, about the microscopic architecture throughout a sample, instead of a local description obtainable with a biopsy. Early demonstrations of microstructure imaging primarily focus on a single imaging device (Alexander et al 2010). However, combining different modalities potentially enhances the technique because different modalities are sensitive to different physical contrasts that may be complementary to one another (Kramer et al 2010).
Disease states frequently alter the microstructure of tissue. An example is the variation of the axon lumen due to demyelination in multiple sclerosis (Lucchinetti et al 2000), or the microstructure degeneration occurring in skin burns (Park et al 2001). Overall, cancer is the most widely known case of disease where different grades (Khan et al 2003) are characterised by different packing of cells, leading to different volume occupied (Enfield and Gibson 2012). For example, nuclei grow with the disease, and occupy larger portions of the cell, or cells might be polynucleate (Shi and King 2005). Different volume fractions occupied by cells, as well as the ratio between nuclei and cell size, are important histological biomarkers (Wittschieber et al 2010). Although histology relies primary on the detection of genetic enzymes, microstructural investigation (using hematoxylin and eosin stain) can also provide an indication of the grade of tumours (Bloom and Richardson 1957, Edge and Compton 2010, Wittschieber et al 2010. Consequently, there is interest for non-invasive histology to obtain clinically equivalent information relaying only on intrinsic contrast (Enfield and Gibson 2012).
One technique used for studying microstructure is diffuse optical spectroscopy (DOS). The microstructure modifies the optical properties (scatter and absorption) of a material and the signal obtained changes accordingly; consequently, it is possible to observe the chemical composition of a sample or the pattern of the microstructure non-invasively (Gibson and Dehghani 2009). The optical absorption (which depends on the different molecules), and the scattering (which depends on the configuration of such molecules or particles) introduces a contrast that can be used to inform a model. DOS provides information about thin emulsions (van Staveren et al 1991, Chantrapornchai et al 1999, Michels et al 2008, food science (Lee et al 1997, Pillonel et al 2003, Everett and Auty 2008, and in astronomy and meteorology (Henyey andGreenstein 1940, Kaplan et al 1994). Other applications of DOS can detect the effects of pathological physiology or microstructure of the tissue, such as change of blood volume, of cytological expression, or of chromophore concentration (Gurjar et al 2001). The sizing of microstructural elements is more commonly implemented with polarized light where cytostructural changes in epithelial cells are observed (Schmitt andKumar 1998, Gurjar et al 2001), but non-polarized light has also been used to distinguish between different microstructures (Ramachandran et al 2007).
A second technique able to investigate the microstructure of a sample is magnetic resonance imaging (MRI), which describes certain chemical and physical properties of tissues. Diffusion MRI provides contrast by the diffusion of water molecules containing spins (Stejskal andTanner 1964, Callaghan andXia 1991). In diffusion MRI, the mobility of spins depends on the density of barriers, the size and permeability of pores, and the intrinsic diffusivity of the medium. Thus, diffusion MRI provides sensitivity to these features of internal microstructure. Diffusion MRI is adopted in food sciences and other manufacturing applications, e.g. in the detection of different properties of emulsions, curds and cheeses (Callaghan and Jolley 1983, Rosenberg et al 1992, Mahdjoub et al 2003, Hollingsworth and Johns 2003, Denkova et al 2004, Johns 2009). It also has major applications in medicine where it is used to detect and diagnose pathologies where histological changes occur, for example in cancer lesions (Koh and Collins 2007), stroke (Neumann-Haefelin et al 1999), multiple sclerosis (Werring et al 1999), and Alzheimer's disease (Rose et al 2000). Recently, researchers have designed biophysical models to relate the signal to microstructural parameters of interest such as a dimension or a permeability of cellular compartment (Assaf et al 2008, Xu et al 2009, Alexander et al 2010. The concept of combining imaging modalities evolved in the last decade, from the basic use of co-registered images down to the estimation of parameters with information obtained from the two modalities. A first approach to multimodality imaging consisted of merging the reconstructed images from each single technique to produce a single output, as Ntziachristos et al (2000) did with the combination of contrast-enhanced MRI and diffuse optical tomography where the simultaneous observation of similar features was used to improve the contrast. A completely different approach was proposed by Huppert et al (2008) who used the raw signals from diffuse optical imaging and functional MRI to improve the estimation of a map of oxyand deoxyhemoglobin distribution in the brain. An elegant approach using the signal from the two modalities to inform a model in order to estimate common information was presented by Atuegwu et al (2010) to estimate the growth of tumors.
Here, we present the first approach based on informing a model with the signal from two complementary modalities in order to obtain a quantitative description of the microstructure. This work demonstrates that microstructural parameters can be estimated with increased accuracy by informing a biophysical model with the signals from two distinct modalities: measurements of the near infrared diffuse optical signal (DOS) and diffusion nuclear magnetic resonance (dNMR). Both techniques can measure the average size of the microstructural pores, but only dNMR can study the distribution of sizes, and only DOS is sensitive to the volume fraction between microstructural compartments where different chemical species compose the compartments. dNMR can provide information about the volume fraction from the peaks of different compounds, but the present implementation is optimized for the study of the oil peak only. This work presents a combined biophysical model able to cross-correct the estimation of the average size by exploiting the contrast mechanism from each modality. The model also improves the accuracy of other parameter estimates as a consequence of a more precise estimation of particle average size.
The article is divided in five sections. After the introduction of the problem (section 1), theory and previous work relevant for this work are presented in the background section (section 2). In section 3, experimental procedures and analytical methods developed are described together with the choice of the specific sample. The results from numerical simulations and experiments are reported and explained in section 4. Finally, a discussion on the different aspects of this multimodality approach is developed in section 5 where advantages and limits of the implementation are discussed.

Diffuse optical spectroscopy
Photon migration in biological media is described by the radiative transfer equation (RTE). The solution of RTE can be obtained with a Monte Carlo approach (Wang et al 1995), the finite element method (Arridge 1995b, Gao and Zhao 2010), or, in a few specific cases, with analytical models. Alternatively, the RTE can be approximated by its polynomial decomposition (P n ), and the solution of first order, the diffusion approximation (DA), can model photon propagation in highly scattering and low absorbing media (Wright et al 2006, Ishimaru 1978. The DA has been developed to represent the propagation of light in the frequency and time-domain (Fantini et al 1994, Contini et al 1997. Here we consider the DA for time-domain, where it has analytical solutions for certain geometries, and one of which describes a transilluminated infinite slab, where the signal of interest is the temporal point spread function (TPSF). An analytical expression of the TPSF, T (t, t 0 ; μ s , μ a , n), is given by Patterson et al (1989) and improved by Patterson et al (1989) and Contini et al (1997). The analytical solution for the time, t, where t is the zero-point of the machine after each pulse, depends on the refractive index, n, on the apparent scattering, μ s , and absorption, μ a , coefficients, on the offset time, t 0 , used for calibration of the experimental signal, and on the time of observation (t) (Contini et al 1997) as with K = 1/(3(μ s + μ a )), v is the speed of light in the medium, z 0 = 1/μ s , z e = 2K, s is the slab thickness, and t 0 is the offset time used for calibration of the experimental signal. Equation (1) is valid for a pencil-shaped source and detector, and this model has been extended (Contini et al 1997) so that a finite detector can be modelled instead of a point detector. The refractive index is fixed a priori, and only the absorption and scattering coefficients are fitted. Mie theory (MT) can relate the macroscopic optical coefficients μ s , and the anisotropy coefficient, g, to the microstructure of the sample. MT is an exact solution of Maxwell's equations for the electromagnetic field around an isolated spherical bead with radius, a (Ishimaru 1978, Wiscombe 1980. MT provides the scattering coefficient μ s,event and the anisotropy coefficient g event , related to a single scattering event as a function of the ratio of internal and dispersion medium refractive indexes (n i and n m respectively), and the shape factor, x = a/(n m λ), where λ is the wavelength of the far field wave. If we consider all the scatterers in the sample, the scattering coefficient of the whole sample becomes where ψ is the volume fraction of the spheres, and the apparent scattering coefficient is μ s = μ s (1 − g event ) (van Staveren et al 1991).
MT assumes that the effect of each sphere on the electromagnetic field is independent of each other. When the spheres are close to each other, MT requires a correction consisting of re-evaluating an apparent volume fraction ψ app (a) which depends on the size of the sphere where ψ (a) is the volume fraction corresponding to size, a, and p the dimensional index, which is 3 in the case of spheres (Twersky 1978). The apparent volume fraction of microspheres with the radius a, when a single radius is present, becomes ψ app corresponding to ψ app,a in equation (3). The other two macroscopic optical properties of the sample, n and μ a , can be related to microstructural properties by assuming they are a linear combination of the properties of the internal and dispersion medium compartments as follows: where μ a,i and μ a,m are the absorption coefficients of the internal compartment and of dispersion medium.

Diffusion NMR theory
Molecules naturally diffuse in a fluid if left free to do so. The dispersion of molecules when no confinements are imposed can be described with Einstein's diffusion equation as the root mean squared displacement where t is the diffusion time, and D the diffusion coefficient (Einstein 1905). The standard diffusion NMR Pulse Gradient Spin Echo (PGSE) sequence applies two gradients of magnitude G and duration δ either side of a diffusion time , interposed by a π refocusing spin (Stejskal and Tanner 1964). For free diffusion, the signal is where b = (γ δG) 2 ( − δ/3), γ is the gyromagnetic ratio, and S 0 is the signal with G = 0 so that an estimate of D can be obtained from two measurements, e.g. one of S and one of S 0 . Identical equations describe a Stimulated Echo (STE).
The microstructure of the sample can introduce barriers that confine diffusion of particles into pores. For example, the dNMR signal from diffusing spins constrained by a closed pore with impermeable walls is approximately (Murday and Cotts 1968) where the coefficients λ m and B m are the structure dependent eigensystem of the Laplacian operator on the constraining geometry. For spherical pores, for example, where a is the radius, μ m is the nth root of μJ 3/2 (μ) − 1 2 J 3/2 (μ) = 0, and J is the Bessel function of the first order. Equations (7) and (8) describe the signal from spherical pores with a single radius; for a distribution of radii the signal is where p is the probability density function of pore radii. Here, following Hollingsworth et al (2003) and Denkova et al (2004), we assume a log-normal probability distribution function for radii with average r and standard deviation σ .

The problem
Here, we consider a simple material in which the microstructure is composed of two compartments: a homogeneous dispersion medium and spherical impermeable pores with the spread of distribution, σ , of sizes, as shown in figure 1. The objective is to evaluate the mean size, r, of the particles in the microstructure, the spread of their size distribution σ , and their volume fraction, ψ (0 ψ 1), using DOS, dNMR, and the combined model.
3.1.1. The sample. The sample adopted in this study is commercial light mayonnaise (Sainsbury's Basic) produced by J Sainsbury plc. Mayonnaise is an emulsion composed of two compartments: rapeseed oil microspheres suspended in a dispersion medium of water and proteins. The microstructure is stable in time, and both DOS and dNMR are sensitive to it. It is isotropic, homogeneous, easy to store and its DOS and dNMR contrasts are similar to those of cells (for dMRI) and nuclei (for DOS) in biological tissues. The components of the emulsion are listed on the label, and they are mainly water, salt (1.5%), egg yolk (6%), and rapeseed oil (27%). Despite the apparent complexity of the dispersion medium, its homogeneity makes it easily usable. The sample has a simple tissue-like microstructure. The distribution of microspheres at the micrometre scales could mimic the distribution of mitocondria and nuclei, with a larger refractive index than the surrounding elements (Subramanian et al 2008). In addition, the fact that the dispersion medium contains proteins and other molecules makes the sample more tissue equivalent than a simple emulsion such as intralipids. From the dNMR point of view, the diffusivity is lower in oil than for water in tissue, but the larger dimension of cells against droplets compensate the differences signals (Xu et al 2009, Panagiotaki et al 2014. However, translation of the technique to biological tissue will introduce various Figure 2. A scheme explaining the algorithm adopted in the combined model approach. From the average size, r, the spread of the distribution of sizes, σ , and the volume fraction, ψ, the DOS and dNMR models generate the synthetic signals, T and S. To define T , Mie theory generates the reduced scattering coefficient, μ s , the absorption coefficient, μ a and the refractive index, n, and T is defined for every time t, with a time offset of t 0 . To obtain S, the diffusion time, , the gradient duration, δ, and the gradient intensity, G are required, and the diffusivity, D is estimated. different effects, such as multiple signal-generating compartments and exchange between those compartments. Nevertheless, simple models along the lines of those we explore here, at least for the dMRI signal, capture the signal fairly well in tumours, see for example (Panagiotaki et al 2014).
3.1.2. The model. The optical model considers a set of homogeneous spherical scatterers with different sizes immersed in a medium. The refractive indices of scatterers and medium are different, and the medium and scatterers (droplets) have different absorption coefficients. The dMRI model considers a single compartment representing the set of spheres inside which spins are confined (droplets) while the external spins are not modelled and they diffuse freely but cannot enter the spherical pores. The parameters of the optical model are the absorption and apparent scattering coefficients of the two compartments, the mean spherical pore radius, and the pore volume fraction. The parameters of the dNMR model are the diffusion coefficients in each compartment and the mean and standard deviation of the pore radius distribution. Thus the models have one parameter in common, r, and several unique for each modality.

Signal models
3.2.1. Diffuse optical signal. We use single-wavelength time-domain DOS measurements to estimate particle sizing, although in literature, it is most common to use frequency-domain measurements. Previous sizing applications use a multispectral approach by analysing a TPSF at a single wavelength to provide sensitivity to μ s and μ a . The limitation to a single wavelength is compensated by the use of a very simple geometry (Michels et al 2008).
The DOS model requires some of the optical properties describing the two compartments as prior information in order to estimate the microstructural parameters. Four parameters are fixed: the absorption coefficients μ a,i and μ a,m , and the refractive indexes of the compartments n i and n m . The oil compartment parameters μ a,i and n i are well known in the literature (van Staveren et al 1991, Michels et al 2008, however the dispersion medium is relatively complex, and changes in the optical parameters n m and μ a,m lead to different estimates of parameters r and ψ. As a first approximation for this work, the parameters selected for λ = 780 nm are reported in table 1 where the dispersion medium is a mix of water (n = 1.33), vinegar (n = 1.31), flour (n = 1.434) and egg yolk (n = 1.419) (Vadehra and Nath 1973). A volumetric weighted average of the different compounds provides the parameters stated in table 1. Similarly for the absorption coefficient, data for the single components were obtained from the literature and an intermediate μ a,m is obtained (Pickering 1992, Chantrapornchai et al 1999. In DOS, the regime is modelled as a diffusive medium, and the microspheres are assumed to be all the same size; therefore, the unique size a is also the mean size of radii r. The signal is approximated by equation (1) in order to obtain the apparent scattering and absorption coefficients. MT is then used to correlate the microstructure with the scattering optical properties, while equation (4) takes refractive index and absorption coefficient into account. Finally, the function implementing the described steps, M DOS is defined. It relates to the microstructure with unique radius r, and volume fraction ψ with its optical properties (μ s , μ a , n) according to MT and equation (4) (μ s , μ a , n) = M DOS (r, ψ ).
where MT provides a relation between r and ψ, and μ s and g, which is corrected with equations (2) and (3), and μ a is obtained from equation (4).

Diffusion NMR.
The sample is modelled as a population of spherical pores (droplets of oil) with a log-normal distribution of sizes, corresponding to the distribution p(a) in equation (9) with mean r, and standard deviation σ . Signals obtained from a STE (Tanner 1970) sequence at different diffusion times and gradients strengths are fitted with equations (7) and (9).
3.2.3. Combined model. The combined model integrates the two previous approaches to obtain a common set of microstructural parameters. Figure 2 presents how the model unifies the information in a common physical model that includes information from both modalities. DOS can provide an estimation of the average radius of microspheres r, and of the volume fraction ψ, while dNMR provides r and the spread of size distribution σ , of the microspheres. In the integration, the DOS model only provides the average radius while dNMR fits the lognormal distribution. Since DOS considers only a single radius of microspheres, the combined model redefines the optical coefficients to model a distribution of radii.
In order to combine DOS and dNMR models, the optical parameters are based on the log-normal probability distribution of sizes, and they are obtained as and where N 0 is the number of microspheres per volume unity, C sca = π r 2 Q sca (a), and Q sca is the ratio between the scattering and the geometrical cross section (van Staveren et al 1991, Michels et al 2008. The final value of μ s = μ s (1 − g), and p DOS (a; r, σ ) is to a log-normal probability distribution function (pdf) of sizes corresponding to the pdf p(a) in section 2.2 (Denkova et al 2004). Note that the DOS data alone do not support estimation of both r and σ so we must represent the distribution of droplets radii by a single mean r parameter. Finally, similarly to equation (10) the optical parameters become (μ s , μ a , n) = M Com (r, σ, ψ ).
3.3. Data acquisition 3.3.1. Diffuse optical signal. TPSFs were recorded in a transmission geometry. A timedomain system, MONSTIR I, was used to generate a non-collimated laser beam of 3 ps pulses at 780 nm with a pulse interval of 12.5 ns (Schmidt et al 2000, Jennions et al 2006. The intensity of the laser beam was reduced with neutral density filters. The mayonnaise sample was poured into a transparent plastic container of 17 × 50 × 50 mm 2 with a wall thickness of 1.5 mm. The laser source was positioned in contact with the external surface of the container, and a detector was aligned on the opposite side, 17 mm away, in order to have a direct transmission measurement. A reference measurement was taken, averaged over 10 s without the sample in position. The entire acquisition required about 1 min. Measurements with the sample in position were acquired by changing the position of source and detector along the sample. All the measures were normalized as T (t, t 0 ; μ s , μ a , n) l 2 = 1. The reference measurement convolved with the modelled signal was iteratively fitted to the measured TPSF. The convolution aims to remove experimental artefacts due to laser drift, fibre coupling and temporal delays in the fibres and cables (Hebden et al 2003).

Diffusion NMR.
We applied a dNMR STE sequence to a 20 × 20 × 20 mm 3 sample using a 9.4T Varian preclinical scanner. A set of measurements with each combination of diffusion times = 50, 100, 150, 200, 250, 400, 700 ms and gradients G = 0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 0.95 T m −1 , δ = 5 ms, and T R = 4 s was obtained. For every set of parameters ( ,δ,G) four repetitions along the three orthogonal directions were taken, and the no diffusion gradient (G 0 ) scan was acquired for each repetition. We acquire the full spectrum in each acquisition and identify only the oil peak corresponding to the signal from inside the spherical pores, which informs the key parameters r and σ .

Parameter estimation
We calculated the microstructural parameters by fitting the model to the experimental data using a Levenberg-Marquardt algorithm. The DOS model estimates r, ψ and t 0 , the dNMR model r, σ and D, while the combined model estimates all the parameters together: r, σ , D, ψ, and t 0 .
The minimization in DOS is performed with a multi-step procedure to enhance the robustness and to avoid convergences to local minima. The parameters are fitted in three steps consisting of different minimizations with Levenberg-Marquardt algorithms initialized with different guesses. The first step uses an arbitrary starting point with physically realistic values for each parameter, the second uses the initial guess of the radius, and the value of concentration obtained by the first step; finally, the third step uses the minimized ψ (from the second step), and the radius fitted in the first step to find the remaining parameters.
The combined model uses the two signals together to define the distribution p DOS (x; r, σ ) = p dNMR (x; r, σ ) = p(x; r, σ ) and the volume fraction ψ. The objective function is a weighted combination of the errors from the two signals (r, σ, ψ ) = min (r,σ ,D,ψ,t where w is the weighting coefficient combining the errors from the two signals in the objective function. A Monte Carlo Markov chain (MCMC) algorithm was adopted to further investigate the properties of the model by sampling the posterior distributions on each parameter after locating the maximum likelihood estimate. MCMC consists of the generation of a Markov chain where the parameters are updated with a Gibbs sampler and the Metropolis algorithm (Walsh 2004, Haario et al 2006. The computations were conducted on a laptop (Quad Core 2.2 GHz, 8GB RAM).

Validation with confocal laser scanning microscopy
To validate the reconstructed microstructural parameters against results describing the microstructure obtained independently, the concentrations of the different compounds were taken from the label on the mayonnaise package, and the size distribution of the oil microspheres were measured using confocal laser scanning microscopy (CLSM).
A small portion of sample was diluted 1:1 with water and stained with Rhodamine 6G, produced by Sigma-Aldrich Co. LLC. This dye stains the proteins surrounding oil globules, and therefore the oil, as red (Blonk and van Aalsh 1993). The diluted mayonnaise was spread on a slide and illuminated with light at 680 nm. A first set of 10 images were acquired varying the position of the focal plane and 5 other images by moving the position of the lenses on the slide in order to capture any differences between areas of the sample. The images were manually processed to determine the diameter of about 100 microspheres. Any spheres visible in the central part of four images were considered, and their diameter measured manually.

Parameters of the models
We first assess the sensitivity of the microstructural parameter estimates to the fixed parameters. Here, the effect of changes in the fixed values of the absorption coefficients μ a,i and μ a,m and the refractive indexes n i and n m , and the weighting factor w are described. Figure 3 shows the change in the estimates with different initial estimates of absorption coefficients and refractive indexes. Synthetic datasets were generated by varying one parameter while keeping the other fixed as in table 1. The offset time t 0 is practically independent of the  (   Figure 4. Effect of the weighting parameter w on the error of fitting, and on the use of synthetic data to obtain the fitting at different w. Here, the parameters considered are the average radius, r, the spread of radii distribution, σ , the volume fraction, ψ, the diffusivity, D, and the time offset, t 0 . In (a) the estimated parameters in a synthetic dataset generated with a SNR = 16, r = 1.5 μm, σ = 1.0 μm, and ψ = 0.3. In (b) the error obtained in the estimation. fixed parameters. The volume fraction ψ remains constant for n m < 1.41 but increases steadily with μ a,m . Similarly, r increases with μ a,m , but a more nonlinear dependence is observed with n m . Figure 4(a) reports the parameters fitted on a noisy synthetic dataset (SNR = 16, r = 1.5 μm, σ = 1.0 μm, and ψ = 0.3) at different values of w. A divergence of the estimates from their real values is observed for w greater than 0.9. Figure 4(b) shows the errors in the estimates at different values of w. Despite the increase in the error as the estimate of r with w, the estimates obtained from DOS are a minimum when w is set to 0.5. The error observed in all the other parameters (also related to dNMR) tends to increase with the value of w. For the remainder of this work, w is fixed at 0.5. Table 2. Estimation of mean radius r in groups of 50 synthetic datasets generated with SNR fixed at 16 and different microstructural parameters with DOS, dNMR, and the combined modality (Com). The datasets are fitted with the two modalities independently and with the combined approach. Here, the estimation of the mean radius r, with the three modalities is applied to datasets generated with r = 1, 2, 3 μm, spread of radii distribution σ = 1 μm, and volume fraction ψ = 0.3. Standard deviations of the parameter estimates are reported in brackets.

Simulation results
The models were first characterized using different synthetic datasets. Tables 2 and 3 show the mean and the associated error (as standard deviation) of the estimates of r and ψ over different datasets. Noise-free synthetic datasets were generated for various combinations of microstructural parameters via Monte Carlo simulations, which are more accurate than the approximate models we fit but computationally more expensive. Noise was added such that the SNR was 16, and the parameters considered were r = 1, 2, 3 μm, σ = 1, 2, 3 μm and ψ = 0.1, 0.2, 0.3. The Monte Carlo code for DOS was implemented by modifying the widely available code MCML (Wang et al 1995) where the tracking of the distance travelled by every single packet of photons in the sample is recorded. The noise is proportional to the square root of the number of events simulated and is Gaussian (Arridge 1995a). The dNMR signal was generated with Monte Carlo simulation (Hall and Alexander 2009) in the Camino toolkit (Cook et al 2006). Table 2 presents the estimates of r. The combined model produces more accurate estimates of r with lower associated error, and the improvement against the single modalities increases with the size of the microspheres. DOS tends to provide better estimations than dNMR in terms of fitted values and error of estimation. The results from the estimation of the parameter ψ are reported in table 3, and an improvement of the estimation with the combined model is observed specifically for lower values of ψ. The identification of σ , directly depends only on dNMR, but is improved by the indirect effect of the enhanced estimation of r by the DOS in the combined approach 4 (see table 4).

Experimental results
The fitting of experimental data is shown in figure 5. Figure 5(a) shows a good fitting of the TPSF since there is a high number of scattering events to ensure the applicability of the DA.  Table 4. Estimation of the spread of radii distribution, σ , in groups of 50 synthetic datasets generated with SNR fixed at 16 and different microstructural parameters. Here, the estimation of σ , with dNMR and combined modality (Com), is applied to datasets generated with σ = 1, 2, 3 μm, volume fraction ψ = 0.3, and mean radius r = 1.5 μm. Standard deviations of the parameter estimates are reported in brackets. Minor artefacts probably due to the reflection in the plastic box occur and are not considered in the model. The dNMR model fits the data well as shown in 5(b) at the longer diffusion times. A lack of fitting precision is obtained for the lowest gradient strengths and diffusion times. The fitted value of the diffusivity coefficient is D = 0.0098(μm) 2 ms −1 , which is in agreement with the expected value found in literature for the diffusivity of fat at room temperature (D = 0.01(μm) 2 ms −1 ) (Denkova et al 2004, Hollingsworth et al 2004.
The estimates of microstructural parameters obtained with the models are reported in table 5 along with the CLSM ground truth. All three parameter estimates from the combined model are closer to the ground truth than from the individual models. DOS appears to overestimate r, while dNMR and the combined model provide a close approximation to the distribution. Although the real value of σ is out of the confidence interval around the estimate from dNMR, an improvement is achieved with the combined model. A slight improvement of the fitting of the volume fraction occurs, but the ability of the method to converge to a suitable final result is highly improved in the combined model. Figure 6 reports the distribution of sizes estimated with the three models compared to the ground truth. It shows a change in the size distribution consistent with the better estimation obtained with the combined model. Figure 7 shows the solutions from a MCMC chain of 10 000 steps after a burn in phase of 1000 iterations for the microstructural parameters of interest and the CLSM results as ground truth. The variance between the elements in the Markov chain indicates the precision of the parameter estimates. The combined model improves the precision compared to the single dNMR model but the results differ from CLSM estimation in terms of σ . At the same time, DOS provides only the mean value of size, which is biased and higher than CLSM. Figure 8 shows the trend of MCMC chains of convergence to the final gradient descent value of parameters. All the Markov chains start from a solution obtained with a single minimization of the objective functions. The convergence is enhanced by the combination of modalities, while the DOS modality requires a longer number of iteration to reach a stable set of parameters. Moreover, the variance obtained with the combined method is reduced compared to the other two models.

Overview of the results
The whole set of parameters may only be obtained by the combination of DOS and dNMR data. DOS estimates r and ψ, and not σ , while dNMR provides r and σ but not ψ. The DOS model requires the absorption coefficients and refractive indexes of the two compartments, which have been retrieved from the literature (Vadehra andNath 1973, van Staveren et al 1991, Figure 7. Monte Carlo Markov Chain solutions in a 10 000 solution chain for the estimation of mean radius, r, volume fraction, ψ, and spread of radii distribution σ .  Pickering 1992, Chantrapornchai et al 1999, Michels et al 2008. The parameters of the inner compartment are well known while the absorption of the dispersion medium has been deduced from its chemical composition. The absorption coefficients and refractive indexes could be obtained more accurately by further experiments, but this problem is outside the scope of this work. Modelling the optical signal with MT introduces an oscillatory dependence of the optical parameters on r (Ishimaru 1978). The objective function of the model introduces multiple local minima identified by similar pairs (μ s , μ a ) but corresponding to different microstructural parameters. The accurate estimation of microstructural parameters depends on a precise choice of the initial guess, and this dependence on initial guess is one of the main disadvantages of this modality. A further disadvantage is the absence of sensitivity to the distribution of sizes resulting from the introduction of the DA, which does not model angular dependence of the scattering.
The dNMR model does not require input parameters other than assuming the shapes and distributions, and the obtained value of the diffusivity coefficient is close to the expected value from literature, confirming the accuracy of the model (Denkova et al 2004). The model is very sensitive to the microstructure and the fitting depends only slightly on the initial guess. The time required for the fitting is around 5-10 s and the model shows a strong sensitivity to the microstructure. However, this modality is not able to provide information about the volume fraction.
The combined model provides enhanced accuracy compared to the single modalities. It is insensitive to the initial guess since it exploits the dNMR properties, but the fitting of distribution is improved thanks to the contribution of DOS, which cross-informs the model about r. Finally, the estimation of concentration is also improved as a result of fitting the size with a robust model such as that of dNMR.

Experimental results
Finding a material able to provide experimental signals from both DOS and was challenging but mayonnaise works well for this simple demonstration. The measurements obtained with DOS were sensitive to the shape of the container, and reflections deteriorated the weak signal produced by the high absorbing sample. Other researchers have used aqueous diffusion medium such as intralipids, but these separate too quickly to be used with dNMR. Another potential confounder to the dNMR acquisitions is the presence of air bubbles.
The results obtained by the use of the combined model on the experimental signals shows an improvement of the fitting of all the microstructural parameters against the single modality models. The comparison against the CLSM indicates a discrepancy between the two measurements due to the spread of the distribution. dNMR is possibly unable to fit the larger microspheres accurately because of to the low diffusivity of fat molecules compared to water (which is more commonly measured with dNMR): even at the largest diffusion time in our experimental protocol, the root mean squared displacement is 3.75 μm. At the same time, a probable coalescence between the smaller bubbles in the sample may produce the narrowing of σ recorded in the CLSM estimation.
The fitting of the parameters related to the size appears to be driven by the dNMR modality. In fact, the distribution of radii estimated with the dNMR and combined modalities do not differ strongly, but the information added in the optical data improves the quality of estimation a little. Since, the optical model is not able to provide an estimation of the spread of the radii distribution, any w factor over-weighting the DOS model will result in instability.

Accuracy of the model
The MCMC algorithm reveals the convergence properties of the three methods. The results presented demonstrate a reduced capability of DOS to define a pair of parameters r and ψ representing the microstructure that generates the experimental TPSF; in fact, many pairs (r, ψ) can provide similar TPSFs and the fitting is not as accurate as with the other modality. MCMC for the DOS modality requires more than 1000 iterations to converge, while dNMR is quicker, and the initial set of parameters estimated by the a single minimization algorithm are very distant from the point of convergence of the chain. The dNMR model is very selective and the parameter estimates are robust to the choice of algorithm. The Markov chains are very regular proving the strong stability of the model, which is able to discriminate between the information contained in different signals.
The MCMC analysis of the combined model illustrates the complementarity of information between the two single modalities. The parameters directly detected by the dNMR model preserve the stability of their chains, showing that properties and advantages of single modalities may be inherited by the combined method. At the same time, the volume fraction estimation is also improved by the interaction between the two modalities. The information contained in the dNMR signal constrains the size estimation and it drives the DOS model to estimate the volume fraction based on an appropriate distribution of radii. At the same time, the mean of distribution of radii is informed by the DOS modality, creating a strong correlation.

Conclusion
The use of combined models increases the accuracy of the fitting because of the interdependence between the single techniques. The approach is general and may be extended to other complementary modalities, e.g. diffraction x-ray and laser scattering signal, by adapting the model only. The demand for a more advanced microscopy investigation, and an improved sample is clear. This approach can be applied to diagnosis of epithelial cancer with the combination of diffusion MR and laser scattering signal.
The two modalities are completely non-invasive. Their contrast is based on intrinsic mechanisms and they do not require contrast agents. Previously, the optical scattering properties have been used to discriminate between different types of microstructure in local measurements (Chung et al 2008) or tomographically, in breast cancer (Wang et al 2005(Wang et al , ,2006). Similar applications of Mie theory on less expensive systems (non time-domain) support the same methodological approach and combination of techniques (Gurjar et al 2001, Oka et al 2011. DOS measurements can be performed quickly and easily, whereas dNMR measurements require a large, static scanner and are slow (around 1.5 h). A combined imaging system may be able to generate good quality fits with a faster imaging system than using dNMR alone. The two modalities are also able to observe different properties of the sample; therefore, they can be used in more advanced models where the correlation between different aspects is important in the description of a pathology.
This first demonstration of the benefits of multimodality for estimating microstructural material features provides motivation for further work on more complex models appropriate for biological tissue and eventually devices enabling multimodal microstructure imaging and non-invasive histology.