Characterisation of heterogeneity and spatial autocorrelation in phase separating mixtures using Moran’s I

In complex colloidal systems, particle-poor regions can develop within particle-rich phases during sedi- mentation or creaming. These particle-poor regions are overlooked by 1D proﬁles, which are typically used to assess particle distributions in a sample. Alternative methods to visualise and quantify these regions are required to better understand phase separation, which is the focus of this paper. Magnetic resonance imaging has been used to monitor the development of compositional heterogeneity in a vesicle-polymer mixture undergoing creaming. T 2 relaxation time maps were used to identify the distribution of vesicles, with vesicle-poor regions exhibiting higher T 2 relaxation times than regions richer in vesicles. Phase separated structures displayed a range of different morphologies and a variety of image analysis methods, including ﬁrst-order statistics, Fourier transformation, grey level co-occurrence matri- ces and Moran’s I spatial autocorrelation, were used to characterise these structures, and quantify their heterogeneity. Of the image analysis techniques used, Moran’s I was found to be the most effective at quantifying the and morphology of phase separation, a quantitative measure comparisons made between a range of phase separation.


Introduction
Many consumer and pharmaceutical products are based on colloidal suspensions [1][2][3], but density differences between the particles and fluid can lead to sedimentation or creaming which will affect the ''shelf-life" of such products. Polymers are frequently added to suspensions to enhance the shelf-life of the product, as well as its efficacy or dispensability [2]. However, adding a nonadsorbing polymer can cause particles to aggregate, via depletion flocculation, causing a gel to form and changing the rate and mechanism of phase separation [3]. By observing suspensions over time, it is possible to better understand phase separation processes, enabling the sedimentation and re-suspendability of a system to be assessed and predicted [1]. For example, creaming in vesiclepolymer mixtures has been studied by using visual inspection [3,4] and optical characterisation [5] to monitor the development of the vesicle-poor, vesicle-rich interface over time. However, visualising phase separation in these suspensions is often difficult, as they are frequently opaque [1].
Magnetic resonance imaging (MRI), however, provides a means by which sedimentation and creaming can be observed, noninvasively, in opaque systems [1,2,[6][7][8][9][10][11][12][13]. Previously, 1 H MR images of spin density have been used to study sedimentation of polymer particles [2,6] and paliperidone palmitate particles [1], and T 2 relaxation time maps have been used to study sedimentation of polymer particles [10,11], glass beads [10] and rayon fibres [9] as well as separation of asphaltenes from crude oil [7] and phase separation in moisturising creams [8]. Chemical shift imaging has also been used to investigate sedimentation and creaming in biodiesel [12] and multinuclear ( 1 H and 19 F) MRI has been used to study sedimentation of polymer particles [13]. In these studies, phase separation has been predominantly probed using vertical profiles of MRI signal intensity, relaxation time or volume fraction, identifying regions either rich or poor in suspended particles. By monitoring these profiles over time, the rate of phase separation can be determined [1,2,7]. However, this approach assumes that the volume fraction of particles only varies vertically, controlled by the direction of gravity, and that there is a complete separation of the different phases, which may not be the case [3,6]. Particlepoor regions within the particle-rich phase have been identified in sedimentation [6] and creaming [3,4] in colloidal suspensions, particularly in suspensions with viscoelastic properties [6]. Yet, these regions are overlooked in 1-dimensional (1D) profiles and analysis, but must be visualised and quantified in order to gain meaningful insight into phase separation processes [6].
By using 2-dimensional (2D) MRI, phase separation can be better monitored and, through image analysis, quantified. An example of such image analysis is segmentation [14], where regions of similar composition can be identified within a sample by comparing the signal intensity or T 2 relaxation time, for a pixel, to a threshold value. This method could provide a means by which the amount of phase separated material may be quantified and can lead to an increase in image contrast and a simplification of features within the MR image [14,15]. Other forms of segmentation include cluster analysis [16], which looks for regions of homogeneous signal intensity, and independent component analysis (ICA) [16,17], which identifies statistically independent groups of pixels. However, the information available by segmentation can be limited and is dependent on the technique chosen. Alternative image analysis techniques include methods based on Fourier transformations or wavelet transformations [18], which determine the spatial frequency of heterogeneity in pixel intensities. These methods can distinguish both fine (high frequency) and coarse (low frequency) features in an image [19]. For example, wavelet transformation allows spatial information to be extracted at specific length scales [20] and can be monitored over time. However, these methods do not easily provide a quantifiable measure by which samples can be compared.
Image analysis methods based on first order statistics have been employed in medical imaging [21], and can be used to describe the distribution of pixel intensities and differentiate between homogeneous and heterogeneous images. Unfortunately, they are unable to provide information about the relative position of features or their connectivity [21]. First order statistics of MRI data have been used to distinguish between healthy and tumorous brain tissue [22]. Other statistical methods include grey level run length matrices (GLRLM) and grey level co-occurrence matrices (GLCM) which assess the probability that specific grey levels occur within a specified spatial relationship and allow the calculation of various parameters, such as local homogeneity, contrast and entropy of the pixel intensities [21,23,24]. However, the number of grey levels chosen in these analyses may affect the result. While using fewer grey levels makes the calculation less computationally demanding, it can also result in a loss of image detail [25].
Image analysis methods using segmentation, transformation and statistical methods have yielded useful information from MR images, however, they can discard useful information about the spatial localisation of the features, be difficult to interpret or over simplify features within an image. An approach used to overcome the limitations of these image analysis methods, is autocorrelation, which quantifies the heterogeneity or clustering in an image. Spatial autocorrelation determines whether an observed variable, at a particular location, is significantly dependant on the value of that variable in a neighbouring region [26]. By quantifying the spatial autocorrelation of a parameter, it is possible to determine whether the data is clustered, as well as quantify how strongly it is clustered [16]. There are a variety of measures of spatial autocorrelation including Moran's I [27], Geary's C [28] and Getis and Ord Gi ⁄ [29], enabling the spatial distribution of a variable to be quantified using a single number. The most widely used of these measures is Moran's I, which has been applied to analyse optical images [30], X-ray CT images [31][32][33] and clinical MR images [16,30,31,34,35]. This method provides a simple means of assessing the degree of spatial autocorrelation, with values ranging from À1, for negative correlations, to 1, for positive correlations. As shown in the images in Fig. 1, increased clustering leads to higher values of Moran's I. Where pixel intensities are randomly distributed, Moran's I is equal to 0, and more alternating features lead to lower, more negative, values of Moran's I [34,36].
In the MRI studies, Moran's I has been used to assess noise levels [30,31,35], study neural networks [16] and distribution of fat in muscles [34]. In the study by Derado et al. [16], Moran's I was used to investigate neural networks, which were identified using segmentation techniques. Spatial clustering of fat, in MR images of muscles, has also been quantified using Moran's I [34]. However, Moran's I has not yet been employed to quantify the amount of heterogeneity within an MR image or characterise the compositional heterogeneity of complex fluids.
In this paper, a vesicle-polymer mixture undergoing creaming was visualised over time using 2D MRI, which revealed vesiclepoor regions with the vesicle-rich phase. The resulting 2D MR images were analysed using first-order statistics, Fourier transformations, GLCM and Moran's I. The results were compared between all image analysis methods and their potential for quantifying phase separation and spatial heterogeneity was assessed. A detailed description of the Moran's I calculation is presented and explained. Different spatial weight matrices were evaluated and it was found that careful selection of the spatial weight matrix made it possible to quantify smaller structures than have been previously accessible [34] using this method.

Theory: Moran's I calculation
Moran's I quantifies spatial clustering and connectivity between clusters by assessing the correlation of parameters between spatially close or contiguous pixels. The more clustered pixels are in an image, the closer the value of Moran's I is to 1 (Fig. 1). Using Eq. (1), Moran's I can be determined for a 2D image containing N pixels, indexed by i and j [36]. Each pixel has a measurable variable, X i/j , and the vector X contains the difference between X i/j and the mean of that measurable variable, X, for the whole image (Eq. (2)). The tensor product, X X, creates a matrix C ij which defines the proximity of i and j, based on Euclidean distance [36]. W ij is a spatial weight matrix, defining the degree of spatial closeness, or contiguity, of pixels i and j [34].
While, the choice of weight matrix influences the value of Moran's I calculated [36], this is frequently not discussed when using Moran's I as a tool for image analysis [16,34]. Yet, the choice of spatial contiguity controls whether the calculated Moran's I value is a measure of global, long-range, short-range or local autocorrelation. Local autocorrelation provides a measure of spatial localisation, which can identify clustering between adjacent pixels within a sample (an image) [37]. Whether pixels are considered adjacent or not can be defined in a number of ways, such as using the Queen's or Rook's case neighbourhood relationships (Fig. 2) [38]. In the case of the Queen's case, correlations between a pixel and its edge and corner sharing pixels are considered, whereas Rook's case only considers correlations between a pixel and the pixels sharing an edge [38]. A spatial weight matrix, W ij , is constructed based on the choice of spatial contiguity (Eq. (3)) and is a binary matrix, where w ij = 1 if pixels i and j are contiguous. The final weight matrix must be symmetric, such that the matrix equals the transpose of itself, W ij = W ij 0 , and w ij = 0 on the diagonal so that a pixel cannot correlate with itself. The cross product of C ij and W ij can be used to determine the final value of Moran's I.

Sample preparation
Dispersions of vesicles formed with the cationic surfactant diethylesterdimethyl ammonium chloride (DEEDMAC) were investigated in this study, which are of particular interest as biocompatible fabric enhancers [39] and have typical vesicle sizes ranging from 100 nm to 5 µm. Aqueous dispersions of DEEDMAC vesicles, with an approximate vesicle volume fraction (U) of 0.5, were prepared according to standard procedures developed by P&G Brussels. Vesicle dispersions were thoroughly mixed with an aqueous solution of the polydiallyldimethyl ammonium chloride (M w = 150000 g mol À1 ) polymer, Merquat 100 (Lubrizol), giving a final polymer concentration of 0.25% wt, and 5 cm tall samples were stored in 20 mm (outer diameter) glass test tubes (Cole-Parmer Instrument Co.). A sample was prepared, which trapped air bubbles during the mixing process. These bubbles were then subsequently extracted using a vacuum pump and then this sample was imaged after t = 21 days. A second sample was prepared in a vacuum, to prevent the trapping of air bubbles, and measurements were made at t = 0 (20 min), 2, 7, 14, 21, 35 and 180 days.

MRI experiments
MRI data were collected on a Bruker DMX 300 spectrometer equipped with a 7 T vertical wide-bore superconducting magnet, operating at a proton resonance frequency of 300.13 MHz with a micro2.5 imaging probe and 25 mm radio frequency (RF) bird cage coil. All experiments were recorded at 293 ± 0.3 K, which was maintained by the temperature of the water cooled micro2.5 gradient coils. The RF pulses were calibrated for each sample. Vertical 1 H MR images of water in the vesicle dispersion were acquired using the spin-echo imaging sequence RARE (Rapid Acquisition with Relaxation Enhancement) [40]. Vertical images were acquired with a 1 mm slice thickness and a 256 Â 128 pixel matrix were used, resulting in square pixels of width 0.156 mm or 0.176 mm. The specific parameters for each image can be found in the I = 0 I I  Supplementary Information (Fig. S1). T 2 relaxation maps were produced from 16 echo images with an echo time of T E = 7 ms and a RARE factor of 16, leading to echo times ranging from 60 to 1740 ms. A repetition time of T R = 10 s was used. All MRI data were analysed using Prospa [41]. A schematic diagram of the slice orientation relative to the sample tube is shown in Fig. 3.

T 2 relaxation time map image analysis
Sections, 50 pixels Â 50 pixels in size, were extracted from T 2 relaxation maps for analysis. These sections were normalised, to account for different mean T 2 relaxation times between the images, by dividing the T 2 relaxation time for each pixel by the mean T 2 relaxation time for the image section. Histogram plots of the population distribution of T 2 relaxation times were generated from each normalised image section. Standard deviation, skewness (c), kurtosis (r) and Sarle's bimodality coefficient (b) (Eq. (4)) [42] were calculated for the normalised image sections using MatLab [43]. 2D Fourier transformations of the selected T 2 relaxation time map sections were generated using Prospa [41]. A grey-level co-occurrence matrix was calculated for each image section, using Phyton, according to Eq. (5), where i and j are the different grey levels and the co-occurrence is dependent on the distance (d) and angle (h) over which co-occurrences are allowed [44]. 256 grey-levels were considered for each image section and co-occurrences were considered between horizontally adjacent pixels (d = 1, h = 90°). Scalar statistics were then calculated for the co-occurrence matrix, which included values for the homogeneity (Eq. (6)), presented here, as well as contrast, correlation and dissimilarity values, which are presented in the Supplementary information [44]. Moran's I values were calculated for each T 2 relaxation time map section, for both the Queen's case and Rook's case neighbourhood relations, using MatLab [43] 4. Results   Fig. 4a shows selected sections from a range of T 2 relaxation time maps for vesicle-polymer mixtures, at different stages of phase separation. The average T 2 relaxation time of water, for the entire section, was in the region 420-490 ms, depending on the sample and its age. T 2 relaxation times were found to distinguish between vesicle-rich or poor regions most effectively, with vesicle-rich regions having values between 370 and 470 ms and vesicle-poor regions between 510 and 800 ms. This is expected to arise from the increase in mobility of water in the vesicle-poor region and decrease in mobility in the vesicle-rich region [45]. The T 2 relaxation time is shorter in vesicle-rich regions, because there is a greater proportion of water molecules inside the vesicles with restricted mobility [45]. The image sections in Fig. 4a(i-iii) and Fig. 4(vi) are of the same sample, which was prepared in a vacuum, to prevent the trapping of air bubbles in the suspension, and were collected immediately after (t = 20 min) the sample was prepared ( Fig. 4a(i)), after 2 days (Fig. 4a(ii)), after 7 days (Fig. 4a(iii)) and after 6 months ( Fig. 4a(vi). For the t = 20 min sample (Fig. 4a (i)), the image appears homogeneous and does not display any evidence of phase separation or the formation of vesicle-poor regions. The image sections in Fig. 4a(ii) and Fig. 4a(iii) were acquired when vesicle-poor regions had formed in the vesicle-rich phase of the sample. The image section in Fig. 4a(vi) was acquired after the sample had fully phase separated into layers and no vesicle-poor regions remained in the vesicle-rich layer. The image sections in Fig. 4a(iv) and Fig. 4(v) are different sections of the same T 2 relaxation time map, which was of a sample which had trapped air bubbles during mixing, that were subsequently extracted using a vacuum pump. The T 2 map was acquired at t = 21 days. The image section in Fig. 4a(iv) displays many, predominantly, vertical vesicle-poor regions in the vesicle-rich phase and the image section in Fig. 4(v) shows the region of the sample where it had phase separated into layers, but with vesicle-poor regions remaining in the vesicle-rich phase. Where the sample developed into two layers ( Fig. 4a(v) and (vi)), the vesicle-poor phase appears beneath the vesicle-rich phase, because the vesicles are less dense than the aqueous continuous phase [4].
Histogram plots, showing the population distribution of T 2 relaxation times for each image section, are show in Fig. 4b. The histogram plots for the first three image sections Fig. 4b(i-iii) are unimodal. The spread of the T 2 relaxation times increases in the plots in Fig. 4b(iv-v), with the plot in Fig. 4b(vi) showing a bimodal distribution. Larger versions of these histograms can be found in the Supplementary information (Fig. S2). Statistical analysis of the data displayed in the histogram plots from Fig. 4b are presented in Fig. 4c-e, giving the standard deviation (Fig. 4c), the skewness and kurtosis values (Fig. 4d) and Sarle's bimodality coefficient (Fig. 4e) for each plot. Fourier transforms of the image sections in Fig. 4a are presented in Fig. 4 f. These are all mostly similar in appearance, except for the Fourier transformation in Fig. 4f(vi), which shows a region of higher intensity running diagonally along the Fourier transformation. The homogeneity values, from the GLCM analysis of image sections (i)-(vi), are shown in Fig. 4 g, with the corresponding contrast, correlation and dissimilarity values in the Supplementary Information (Fig. S3). Finally, in Fig. 4 h, the values for the Moran's I, calculated using both Queen's and Rook's case neighbourhood relations, are presented for each for image section in Fig. 4a.   Fig. 5a shows T 2 relaxation time maps for a single sample over time, giving the Moran's I value for each image. Vesicle-poor regions can be seen to develop in this sample after 2 days, with the sample splitting into two layers after 14 days. However, it is only after 35 days that the vesicle-poor phase can be observed by eye. The value of Moran's I is plotted as a function of time in Fig. 5b and there is a general trend of increasing Moran's I as the sample becomes more phase separated over time. In Fig. 5c, a plot of the natural logarithm of Moran's I against time (Fig. 5c) in presented, which shows a linear relationship at the earlier time points ( 35 days).

Discussion
From the T 2 relaxation time maps in Fig. 3a, it is apparent that not all the systems separate into two distinct layers. Also, of the two examples, where two layers can be observed (Fig. 3a(v and  vi)), only one of them is fully separated, with one of them still showing significant heterogeneity in the top layer (Fig. 3a(v)). Hence, 2D imaging and image analysis is essential to detect and quantify phase separation in these vesicle-polymer mixtures. From the histogram plots (Fig. 4b), showing the T 2 relaxation time distribution for each image section (Fig. 4a), it can be seen that the formation of vesicle-poor regions only becomes apparent once the there is significant phase separation (Fig. 4a(iv-vi)). While, the presence of smaller, partial, phase separation (Fig. 4b(ii) and (iii)) cannot be detected in their corresponding histogram plots, which look almost identical to the one for the homogeneous sample ( Fig. 4b(i)). Similarly, the values calculated for the standard deviation and Sarle's bimodality coefficient also showed little difference between the first three image sections. It is only when larger vesicle-poor regions are formed or when different layers form that this can be detected using the standard deviation or Sarle's bimodality coefficient (image sections (iv-vi)). The skewness and kurtosis values were found to vary the most between systems, but not in a unique manner, preventing a means of quantifying the developing phase separation. For example, the features observed in image sections (ii) and (iv) gave similar kurtosis values, while the image sections (iv) and (vi) had similar skewness values. However, all three images showed very different degrees and morphologies of phase separation. The Fourier transforms were less able to distinguish between the different image sections, unable to identify the change in the spatial frequency of sample heterogeneity and did not offer a simple quantifiable measure of phase separation. From the GLCM analysis of the image sections, only the homogeneity values were found to show any correlation with the extent of phase separation, with there being a general trend of an increasing homogeneity value as the phase separation progressed. While the GLCM homogeneity value was able to distinguish between the slight phase separation observed in image sections (ii) and (iii), from the homogeneous image section in (i), it was unable to distinguish between different degrees of phase separation, such as those observed for image sections (ii) and (iii) from that observed in image section (iv).
It is only the values of Moran's I, shown in Fig. 4 g, that provide a means of quantifying the heterogeneity of phase separation, where the degree of phase separation appears to correlate with the value of Moran's I. Moran's I is lowest for the most homogeneous sample, where no phase separation is detected, but steadily increases as the area of the vesicle-poor regions increases. It is only the Moran's I value that is able to detect the small heterogeneities in image sections (ii) and (iii) and monitor the evolution of phase separation. The Moran's I values for image sections (v) and (vi) are the highest, as all pixels with the higher T 2 values, in the lower phase, correlate strongly with each other and there is little variation in T 2 in this region. It is found that, in certain situations, the value of Moran's I depends on which weight matrix is used, as demonstrated in the case of the two image sections showing the least phase separation ((ii) and (iii)). These weight matrices have enabled smaller features to be detected by Moran's I, than has been previously observed [34], and where the Rook's case weight matrix is used, the values for Moran's I are higher, suggesting it is the most effective at detecting small heterogeneities and distinguishing between samples at the early stages of phase separation. This is demonstrated in the time series of images in Fig. 5a, which show the development of phase separation before the sample has formed the two distinct layers (at 35 days) that can be observed visually. The plot of Moran's I vs time indicates an exponential relationship, which is confirmed in the plot of ln (Moran's I) vs time. As there is a linear relationship between ln (Moran's I) vs time, at early time points ( 35 days), the rate of phase separation appears to follow first-order kinetics. This has been previously observed during the creaming process of other vesicle-polymer mixtures [3]. In these studies on the kinetics of phase separation, the process has been described as poroelastic consolidation [3]. From the plot in Fig. 5c, a rate constant of 0.002 h À1 was determined, which leads to a value for the timescale for collapse, s, of 500 h. This timescale is believed to be dependent on viscosity, particle volume fraction, initial sample height, network permeability and the elastic strength of the gel [3,46]. A previous study of a similar system found a s value of 20-180 h, for samples of comparable height, though at lower vesicle volume fraction and with a smaller polymer, [3]. However, this study was only able to observe the creaming process once the vesicle-poor lower phase had formed, which prevented the kinetics of this process being characterised in the early stages. Thus, we have shown, in this paper, that by combining 2D MRI, with Moran's I analysis, it is possible to detect and quantify phase separation, immediately after onset, and hence gain insight into the kinetics of these processes much earlier than through visual observation.

Conclusion
Sedimentation and creaming of colloidal suspensions can lead to particle-poor regions in the particle-rich phase [3,4,6]. These regions are ignored when sedimentation is considered as a 1D process, controlled only by gravity and hence, 1D imaging and analysis [1,2,7] is not always appropriate. In this paper, we have shown that particle-poor regions can form in the particle-rich phase, during creaming in vesicle-polymer mixtures by using 2D T 2 relaxation time maps. By observing creaming in a range of vesicle-polymer mixtures, over time, and applying a variety of image analysis techniques, it has been found that Moran's I is the most effective at quantifying compositional heterogeneity caused by phase separation. Moran's I was the only technique able to differentiate a homogeneous, pre-phase separation, image from an image displaying small heterogeneities, at the early stages of phase separation and able to quantify the degree of phase separation. This sensitivity to small heterogeneities is especially important when quantifying phase separation which is not characterised by a vertical variation in volume fraction [3,4,6]. As such, Moran's I has a provided a robust and sensitive means of quantifying phase separation, which enables the relationships between sample composition, environmental history and time to be correlated. The importance of the choice of weight matrix was demonstrated as it was found that by choosing the right weight matrix made it possible to quantify very small heterogeneities more effectively. Moran's I was also used to follow phase separation over time and it was found that the phase separation process followed first order kinetics as expected for a colloidal gel undergoing poroelastic consolidation [3,46]. However, by using Moran's I to quantify the heterogeneity in 2D T 2 relaxation maps, information about the rate of phase separation was obtained prior to formation of the lower, vesicle-poor, phase. For all these reasons, Moran's I has been shown to be an invaluable tool in the study, and improved prediction, of highly complex, multi-factor, sedimentation or creaming of colloidal suspensions. To this end, we believe greater insight is possible, into the mechanisms controlling creaming and sedimentation in vesicle-polymer mixtures, using Moran's I.