The application of non-linear curve fitting routines to the analysis of mid-infrared images obtained from single polymeric microparticles

For the first time, we report a series of time resolved images of a single PLGA microparticle undergoing hydrolysis at 70 °C that have been obtained using attenuated total reflectance-Fourier transform infrared spectroscopic (ATR-FTIR) imaging. A novel partially supervised non-linear curve fitting (NLCF) tool was developed to identify and fit peaks to the infrared spectrum obtained from each pixel within the 64 × 64 array. The output from the NLCF was evaluated by comparison with a traditional peak height (PH) data analysis approach and multivariate curve resolution alternating least squares (MCR-ALS) analysis for the same images, in order to understand the limitations and advantages of the NLCF methodology. The NLCF method was shown to facilitate consistent spatial resolution enhancement as defined using the step-edge approach on dry microparticle images when compared to images derived from both PH measurements and MCR-ALS. The NLCF method was shown to improve both the S/N and sharpness of images obtained during an evolving experiment, providing a better insight into the magnitude of hydration layers and particle dimension changes during hydrolysis. The NLCF approach facilitated the calculation of hydrolysis rate constants for both the glycolic (kG) and lactic (kL) acid segments of the PLGA copolymer. This represents a real advantage over MCR-ALS which could not distinguish between the two segments due to colinearity within the data. The NLCF approach made it possible to calculate the hydrolysis rate constants from a single pixel, unlike the peak height data analysis approach which suffered from poor S/N at each pixel. These findings show the potential value of applying NLCF to the study of real-time chemical processes at the micron scale, assisting in the understanding of the mechanisms of chemical processes that occur within microparticles and enhancing the value of the mid-IR ATR analysis.


Introduction
Multivariate analysis of hyperspectral imaging data is a rapidly developing research eld and has received considerable attention over the last decade. 1 Hyperspectral images that provide spatial and spectral information at the same time in an array of pixels and an individual pixel, respectively, can be obtained by a number of techniques including, X-ray tomography, 2 X-ray uorescence, 3 Raman microscopy, 4 and near infrared (NIR) imaging. 5 Over the last two decades Fourier transform infrared (FTIR) spectroscopic imaging has become routine, facilitating chemical characterization of multicomponent systems under both static and kinetic conditions. 6,7 Currently more than 80% of all pharmaceutical formulations are delivered in a powder format 8 and attenuated total reectance (ATR) mode has proved advantageous, particularly for pharmaceutically relevant systems, because the sample can be in any phase, form or shape therefore no sample preparation is necessary. 9 Micro-ATR-FTIR imaging with a germanium objective provides a higher spatial resolution compared to transmission and transection due to the effective magnication imparted by the high refractive index of the ATR crystal material. Conveniently, ATR-FTIR imaging in macro mode, i.e. without the use of a microscope, provides a temperature controlled environment for studying dynamic systems with a larger eld of view. 10 In the macro ATR mode, kinetic processes can be probed with IR light such that each of the 2D array of pixels of the focal plane array detector acts as an individual detector, allowing the collection of thousands of IR spectra simultaneously. Consequently, a stack of 2D images at a range of IR wavelengths can be collected within a few minutes, proving good temporal resolution for relatively slow processes. 11 The use of ATR sampling in infrared spectroscopy is based upon the fact that, although total internal reection occurs at the sample-crystal interface, radiation does in fact penetrate a short distance into the sample, this is known as the evanescent eld. The distance that the evanescent eld can travel within a sample in direct contact with the ATR crystal is dened as the depth of penetration (d p ). Harrick and duPre 12 dened d p as the value at which the initial electric eld strength (E o ) decays to a value of E o exp À1 and can be given as; where l is the wavelength of light, q is the angle of incidence, n 1 and n 2 are the refractive indices of the ATR crystal and the sample in contact with it, respectively. 13 For example, for the conditions of the experimental setup and sample used in this study; diamond crystal (n 1 ¼ 2.42) in contact with PLGA (n 2 ¼ $1.45), eqn (1) gives a calculated d p of $1 mm for l ¼ 5.7 mm (1745 cm À1 where the PLGA carbonyl band shows high absorption). In practice, the true depth of penetration is $3 times more than the calculated d p value. 13 Clearly, ATR imaging only probes the near-surface of a sample but this permits the study of samples in aqueous media (n 2 ¼ $1.33) which can be very challenging using traditional approaches such as transmission. Eqn (1) also shows that a good optical contact between the sample and ATR crystal is critical for obtaining a uniform absorbance prole avoiding artefacts within the eld of view. The Golden Gate™ Imaging Single Reection Diamond ATR Accessory (Specac Ltd) has corrective optics that adjust the plane of best focus to be situated on the crystal surface thus minimising any anamorphism and an autolevelling sapphire anvil that ensures a uniform contact between the crystal and sample. 14 FTIR images, consisting of spectral and spatial information, must be analysed in detail to convert the data into chemically and physically signicant information. 6 Several methods for spectral and spatiotemporal data modelling can be used independently or in combination, however, the choice of analysis approach(es) for a series of hyperspectral image sets is determined by the nature and quality of the spectra and the information that needs to be extracted. 15 The main challenge in hyperspectral image analysis on timeresolved multicomponent data sets is the extraction of the important information from the large volume of data generated including overlapping spectral features and noise. For time resolved experiments that continue for longer than a few hours unavoidable contributions from variations of atmospheric water vapour during the experiment has a huge effect particularly when the spectral region of interest is between 1500-1700 cm À1 . This creates complications in the image analysis because the true peak centre, required for peak height measurements, may vary between pixels due to the superimposition of the rotational ne structure of atmospheric water onto the sample spectrum. Algorithms for subtracting atmospheric water vapour superimposed on the sample spectrum using a known water vapour spectrum can be employed for bulk ATR measurements collected using a single pixel detector. However, this does not always eliminate the problem in hyperspectral imaging, because of spectral variations between pixels emerging from the inherent low signal to noise ratio (S/R) compared to single point IR measurements. Applying a derivative is a common approach which eliminates baseline effects in univariate analysis but the water vapour bands are magnied by derivatives to such a degree that spectral information from the sample may be difficult to observe. Therefore spectra need to be deconvoluted by using so or hard multivariate methods so that bands from the sample can be elucidated from the interfering water vapour signal. The most commonly applied multivariate tools to extract information from hyperspectral images include principal components analysis (PCA) and multivariate curve resolution (MCR) and comparisons between their application have been discussed elsewhere. 1,16,17 The aim of this paper is to evaluate the advantages and limitations of univariate, hard and so multivariate approaches for the analysis of ATR-FTIR images, collected during a dynamic process in real time. We also investigate the suitability of each of these analysis approaches for extracting quantitative information regarding the changes in the hydration and chemistry of a single poly(lactic-co-glycolic) acid (PLGA) microparticle and calculate the reaction rate during hydrolytic degradation.

Experimental
Materials PLGA RG752H (75 : 25 lactide : glycolide, I.V. 0.16-0.24, Bohringer-Ingleheim) pharmaceutical grade CO 2 (BOC Special Gasses) were used as received. Water used in the experiments was puried with the ELGA Purelab option-R water distillation apparatus (Up to 15 MU cm, Type II water) and degassed using a Fisherbrand FB11004 ultrasonic bath at ambient temperature and 100% ultrasound power for 15 minutes.

CriticalMix™ process
PLGA is a random copolymer of poly(glycolic acid) (PGA), poly(lactic acid) (PLA), and is a U.S. Food and Drug Administration (FDA) approved, biodegradable 18,19 synthetic polyester that is physically strong and highly processable. 20 PLGA has suitable properties for biomedical applications as a scaffold 21 and sustained release systems 22 and has been comprehensively studied as carrier matrix for macromolecules such as proteins and peptides which are considered promising for the treatment of a range of conditions such as cancer, human growth deciency, and multiple sclerosis. 23 The manufacturing process used to produce the microparticles used in this study was a simple, one-step process which has been used recently to encapsulate protein based drugs with 100% encapsulation efficiency. 22,24 When PLGA is exposed to scCO 2 in a pressure vessel, the polymer is liqueed. If the liquid is depressurised through a nozzle, whereby the CO 2 returns to a gaseous state, the polymer solidies resulting in the production of microparticles. Careful control of the nozzle dimensions and depressurisation rate determines the particle size. Batches used in this study were prepared by adding 2.1 g of pre-weighed polymer, to the Particles from Gas Saturated Solutions (PGSS) apparatus.
The apparatus was sealed, pressurised with CO 2 to 700 psi (48 bar) and heated to 40 C. Once at temperature, the pressure was increased to 2030 psi (140 bar). The liqueed polymer was then stirred at 150 rpm for 1 hour, aer which time stirring was stopped and the mixture was depressurised through a nozzle generating microparticles. These were collected in a cyclone and recovered as a free owing white powder.
Real time ATR-FTIR imaging of reactions Infrared images were collected using an Agilent 680-IR FT-IR spectrometer coupled with a liquid nitrogen cooled mercury cadmium telluride focal plane array detector MCT-FPA (64 Â 64 pixels), capable of simultaneously collecting 4096 spectra from an image area of 640 mm Â 640 mm using the Golden Gate™ Imaging Single Reection ATR Accessory (Specac Ltd). This accessory has a diamond internal reection element with corrective optics that adjusts the plane of best focus to sit on the crystal surface, eliminating any distortion, so that a symmetrical point spread function can be assumed. The angle of incidence of the infrared beam was 45 and the numerical aperture (NA) of the system was 0.32. Images were recorded in rapid scan mode and typical collection times were $5 minutes. The detector was mounted on a Large Sample (LS) external sample compartment and the infrared beam from the spectrometer was projected directly on to the FPA aer passing through the ATR sampling accessory. The physical size of each FPA pixel is 40 mm Â 40 mm and with 4Â magnication, each pixel represents a 10 mm Â 10 mm square in the image.
To set up the hydrolysis experiment, a single PLGA microparticle was placed in direct contact with the ATR crystal using a 40Â microscope standing on top of the ATR accessory and sufficient pressure was applied using the sapphire anvil to ensure good contact between the particle and the crystal resulting in some deformation of the spherical particle. The images collected were not circular, most likely showing some evidence of anamorphism despite the use of corrective optics and this issue when using such collection optics has been observed previously by Everall et al. 25 and Chan et al. 26 Once a satisfactory 'dry' image was collected water was introduced into the chamber in such a way that access to the particle was limited to the sides only as depicted in Fig. 1. Images were collected using the Agilent Technologies' ResolutionsPro FTIR Spectroscopy soware version 5.2.0(CD846) at pre-determined time intervals and the collection parameters used were 128 co-added scans at a 4 cm À1 spectral resolution, in the mid-infrared (MIR) range (3800 to 950 cm À1 ). The raw processed images were obtained by ratioing against a background of the blank ATR crystal comprising of 256 co-added scans.

Data pre-processing
Raw processed image les were cropped between 1820 and 1000 cm À1 which provided a number of characteristic bands associated with PLGA and also included the water d(OH) peak $1635 cm À1 . In order to remove systematic discrepancies such as variations in detector sensitivity or sample thickness, raw images were vector normalised to minimise absorbance variance. Vector normalisation works such that each spectrum is divided to the vector length which is square root of the sum of all absorbance values squared. Instrumental factors such as uctuations caused by the changes in the IR source intensity or temperature or detector sensitivity and optical artefacts caused by mismatch of the refractive index of species being imaged are known to cause a slope in the baseline in FTIR spectroscopy. In order to eliminate this slope, a second derivative can be applied, however caution is advised if there are very weak infrared bands of interest as a second derivative will magnify the noise and some low intensity peaks may be lost in the noise. However a rst order baseline correction between the two end points of the previously cropped spectral range was found useful for elimination of the baseline dri when collecting an image. The full data pre-processing procedure is described in the ESI section including Fig. S1. † Univariate analysis. Univariate analysis, or functional group imaging, only considers the peak height of, or the integrated absorbance under, a peak of interest, therefore chemical information can be obtained based on the association of the peak position with certain functional groups. Although it is useful in providing a quick overview of the species in the raw image data, a peak height image usually convolutes the underlying chemical information from the overlapping peaks in Fig. 1 (a) A microparticle on the Agar microtool straight needle T5340 to be placed on the ATR crystal. (b) A microparticle placed on the ATR crystal before the anvil was brought into contact with it. (c) Schematic of the experiment setup after the particle was placed, anvil was in contact and water was added.
This journal is © The Royal Society of Chemistry 2014 Analyst, 2014, 139, 2355-2369 | 2357 the intensity map. Integrating the absorbance for a particular peak over a spectral region of interest may provide a better distribution map by increasing the S/N. However, as in most cases, a complex heterogeneous mixture will not include an isolated band of interest therefore the integrated values will not represent the amounts present at different locations in the image with the required certainty. Another draw-back of functional group imaging is that a homogeneous object may be shown to be in a different location in the image when univariate images are created based on the integration of infrared bands observed in different parts of the spectrum. This occurs because the diffraction limited spatial resolution of the generated images strongly depend on the chosen part of the MIR spectral region. 27 These disadvantages can be overcome by using multivariate approaches, which allow all of the species to be searched in the same spectral range thus allowing a contribution from all the peaks of the same species (or factor). Another advantage of multivariate approaches is that the interferences such as water vapour can be detected as a factor, and may therefore be eliminated.
Multivariate curve resolution-alternating least squares. Multivariate curve resolution-alternating least squares (MCR-ALS) is a tool that facilitates the extraction of information from various spectroscopic and imaging techniques and is widely used in physical and biological sciences. 28 MCR methods are so modelling tools that require no prior knowledge of the nature of the components in the mixtures to be deconvoluted 29 and can be grouped into, non-iterative and iterative approaches. Non-iterative MCR algorithms are unconstrained therefore may nd unique proles within the mixture being studied. However although these mathematically obtained unique proles may be close to real chemical proles that are to be resolved from the mixture data set, non-iterative methods oen suffer from ambiguities that arise due to strong overlap between components of the species or low S/N of the data. On the other hand, with iterative MCR algorithms, so concentration constraints such as non-negativity and unimodality or hard constraints such as the use of pure component proles, which are compared with extracted proles and modied accordingly, are implemented to minimise ambiguities. 30 By denition, MCR-ALS is a somodelling method that can extract component information from the raw measurement data alone, as long as this data contains some variance; spatially as one might anticipate in an image or as a function of time when monitoring a reaction. Detailed information about curve resolution techniques can be found elsewhere 31 however we would like to provide the necessary theory on the MCR-ALS soware employed here (MCRv1.6) developed by Andrew and Hancewicz. 32 The MCR-ALS algorithm used in this paper 32 has been developed for two-way data, therefore to use this approach, the 3D hyperspectral data cube must be unfolded into a twodimensional matrix and refolded aer analysis as shown in Fig. 2.
FTIR images have a bilinear structure and it is assumed that some form of a Beer-Lambert relationship exists between spectral intensity and concentration. 33 Therefore a bilinear model representing a spectroscopic image can be given as, where D is the measured data matrix, A is the matrix of normalised spectra of pure chemical components and B is the related intensity matrix for each component. 34 The matrix size indices, u, v and z represent the number of spectroscopic resolution elements (wavenumbers), total number of spectra and number of resolvable components, respectively. Rearranging eqn (2) for the least squares estimation of A and B yields: To generate optimal matrices of A and B, one rst needs to estimate their initial values. In this paper, a non-linear iterative partial least squares (NIPALS) decomposition method has been applied, the merits of which when compared to using random numbers, eigenvalue decomposition, or dissimilarity criterion is discussed in detail elsewhere. 33 This is followed by the selection of the optimal number of factors to calculate and is achieved by considering the appropriateness of the initial estimate of the number of components. The nal step, factor rotation, is a rening process where alternating least squares (ALS) is used to determine the optimal loadings and abstract factors matrices such that once recombined they most closely resemble the hyperspectral data matrix.
Although ALS is the most common method to decompose eqn (2) iteratively, a modied alternating least squares method (MALS) proposed by Wang et al., 34 which has been shown to overcome unstable convergence properties giving a nonoptimum least squares solution, has been used for decomposition for all the MCR images generated here. The MCR-ALSv1.6 soware was run from its graphical user interface (GUI) that allowed the user to input the number of factors (or components) to be estimated, number of iterations (used 500), and an ALS non negativity constraint (in this instance we used MALS-2D). The rationale for the use of 2 factors is justied based on the fact that the image quality (S/N, resolution) is better when using 2 factors than when using 3 or 4 factors and the same as using 5 factors, see ESI, Fig. S2. † And it is important to note that MCR results might be improved by using different strategies that were not covered in this paper and/or more so constraints or hard constraints. 35,36 For example, when the data were cropped down to the 1820-1500 cm À1 region which includes the PLGA carbonyl ($1745 cm À1 ) and water dOH ($1635 cm À1 ) MCR results were improved and although LA and GA peaks were not deconvoluted using MCR even when used cropped data between 1500 and 1300 cm À1 , they were deconvoluted when 2 other PLGA compositions (PLGA100/0 and PLGA50/50) were included in the data set.
Another interesting, but rather lengthy, strategy that could have been considered was removing the solvent IR background contributions as developed and demonstrated by Kuligowski et al. 37 in liquid chromatography infrared detection. They have estimated the solvent background successfully and by two different methods; principal component analysis and simple-to-use interactive self-modeling analysis (SIMPLISMA) and aer subtracting the estimated background contributions from their data sets, MCR-ALS provided improved S/N ratios and resolved overlapping chromatographic peaks. As the shape of the d(OH) band has been shown to change from pixel to pixel in FTIR imaging data, we considered that using the simple subtraction of a pure water spectrum at each pixel would introduce spurious peaks due to imperfect subtraction and hinder the MCR analysis. Therefore, for the purposes of this study and to keep the MCR analysis as a so modelling approach, we decided to use no prior information of the pure components and utilise the whole ngerprint region of the IR spectra (1820-1000 cm À1 ) for the MCR analysis. We have not applied hard constraints such as pure spectra and/or different compositions of pure PLGA spectra as initial estimates of the MCR factors in the image set. Each of these would have been valid approaches but were beyond the scope of this work.
Nonlinear curve-tting. The strategy behind curve tting or hard modelling is to generate a model spectrum based on the sum of the component peaks that it contains. If the components within a spectrum and by extension an image are known, then an estimate of the relative amount of each component can be used as a starting point in the application of hard modelling. Here we have developed an efficient curve-tting algorithm to optimise the parameters (lineshape, peak height, peak width) for infrared absorption bands and the percentage contribution of each tted curve to the overall spectrum was used as the initial loading value.
In order to use curve-tting procedures, analytical functions must be used which describe the lineshapes of the peaks. Typical lineshapes encountered in spectroscopic studies are the Gaussian and the Lorentzian where I(k) is the intensity at wavenumber k, k 0 is the wavenumber at the peak centre and D is the full width at half maximum (FWHM). 38 The Lorentzian peak shape is oen used to t infrared absorption bands. 39 However in real infrared spectra, effects including hydrogen bonding, rotational ne structure and, in FTIR-ATR, anomalous dispersion effects might affect the shape of the infrared band thus a true Lorentzian shape does not always occur. Due to the asymmetry and complexity of the bands within infrared spectra, these standard lineshapes should only be used with caution. In such situations a third lineshape function, the Voigt prole representing the convolution of a Gaussian and a Lorentzian, is oen used. 39 However in practice, the use of the exact Voigt prole in curve tting can be time consuming because it involves repeated convolutions. Each of the mentioned lineshapes are symmetric functions, but it is very important to note that unlike traditional IR spectra, IR imaging spectra are more susceptible to asymmetry even in transmission mode. It has oen been assumed that recorded infrared imaging data are exactly like their traditional single point counterparts of bulk materials, and the spatial geometry of the sample was not thought to be an important factor in the processed IR image. 6 In transmission and reection absorption measurements, it is recognized that there are differences in the data recorded from a bulk measurement and a microscopic measurement, caused by light focusing at the point of interaction with the sample when the microstructure of the sample is of the same length-scale as the wavelength of the interrogating radiation 40 and differences in the refractive indices of the rarer media at air/sample interface can result in artefacts in mid-IR images. In ATR-FTIR dispersion effects can lead to asymmetry in observed infrared bands, therefore in this paper we have considered the use of a family of peak curves called the Pearson proles, specically the Pearson IV prole which is asymmetric and is related to the Pearson VII, where P ¼ ½2ðk À k 0 Þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1=M À 1 p =D and M is known as the Pearson parameter. Eqn (7) reduces to the Lorentzian function when M ¼ 1, approaches the Gaussian function when M becomes large and can approximate to the Voigt function for intermediate M. 38 Modifying eqn (7) with an exponential term provides the required asymmetry and for M < 1 the distribution has very broad wings. The Pearson IV function is given in terms of this by The exponential term has been shown to affect peak shapes when M is close to 1 (i.e. Lorentzian) more than it does for those where M is close to 10 (i.e. Gaussian). 38 A routine in MATLAB Version 7.10 (R2010a), was developed to create ts for the each pixel within an IR imaging data set consisting of 4096 spectra individually. Peaks were detected using the algorithm described elsewhere 41 and height and width of the peaks in the raw spectra at detected positions were used as initial guesses for the iterative loop. Application of non-linear curve-tting to large two-dimensional experimental infrared spectroscopic arrays 39 and a combination of non-linear curve-tting and self-modelling curve resolution (SMCR) 42 have been demonstrated. Soware packages such as PyMCA, which is a non-linear least-squares tting application have previously been developed for X-ray imaging data. 43 However to the best of our knowledge, we are demonstrating a peak detecting non-linear optimisation algorithm applied to temporal ATR-FTIR imaging data for the rst time.
The method by which we have approached this tting procedure is described in detail in the ESI † and examples of the output from of the peak detection algorithm are given in Fig. S3. † To summarise our approach: 1. We select 18 spectra from the image based on their position; 6 from the water rich regions, 6 from the PLGA rich regions and 6 from the interface.
2. We iteratively and in a user supervised manner, adjust a peak detection threshold, such that all peaks in all 18 spectra are automatically detected by the peak detection algorithm. We x the peak detection sensitivity parameter 'Delta'.
3. Using the xed peak detection parameter we apply the peak detection algorithm to all 4096 spectra in an unsupervised manner.
4. We then t peaks of all 4096 spectra based on the number of detected peaks and their positions using a gradient search algorithm which nds the best t for the measured data optimising the variables; peak height, peak width, peak centre, Pearson parameter and asymmetry parameter for each peak. 5. We then use the intensity of the peak $1635 cm À1 (selected in a supervised manner) to determine the distribution of water in the system and the sum of the intensity of all the peaks except the peak $1635 cm À1 to determine the distribution of PLGA.
6. The intensity of the peaks $1456 and $1424 cm À1 , were used to calculate k values for lactic acid and glycolic acid units respectively. Fig. 3 shows the result of processing the same infrared image of a single PLGA microparticle in four different ways; using the peak height of a single peak in this case the ester carbonyl at $1745 cm À1 (Fig. 3(a)), plotting the distribution of a factor identied as deriving from PLGA in MCR-ALS (Fig. 3(b)), by plotting the sum of 10 peaks tted between 1800 and 1000 cm À1 which does not include the peak $1635 cm À1 assigned to the water bending mode resulting from a NLCF procedure (Fig. 3(c)) and by plotting the sum of 10 peaks tted between 1800 and 1000 cm À1 which does not include the peak $1635 cm À1 resulting from a peak tting procedure using solely Gaussian lineshapes procedure (Fig. 3(d)). Images were generated from measurements conducted on an as received PLGA microparticle and on the same PLGA microparticle, immediately aer water had been brought into contact with it.

Spatial resolution comparison
The spatial resolution of a microscope is theoretically determined by the diffraction of radiation i.e. the Rayleigh criterion, which is dened in eqn (9) as where l is the wavelength and NA is the numerical aperture. This implies that two objects are totally resolved if they are separated by 2r. 44 Under these conditions, 2r is the denition of spatial resolution. However in an infrared imaging system, FPA detector pixels are not points and have a nite size that is greater than the wavelength of the IR light, therefore this relation is never observed. Furthermore, in the ATR experiment, the penetration depth (eqn (1)) can degrade the lateral resolution. Therefore for FPA imaging systems it has been shown to be more appropriate to determine the spatial resolution based on real measurements. 25,26,45 The practical method we have used is the 'step-edge' method 25 which is based on the observation of a step change increase when the intensity of a selected wavelength is plotted along a chosen line parallel to one of the axes in the 2D image. This step shape represents the Line Spread Function (LSF) (Fig. 3). The derivative of the LSF with respect to its variable (which is position or pixel number) is described as the Point Spread Function (PSF). The PSF is the response of the system to a point source and is generally considered to be an Airy function. The FWHM of the PSF gives the spatial resolution of the imaging system. Therefore, the FWHM of a Gaussian that is tted to the PSF has been dened as the spatial resolution in this paper as described by Offroy et al. 44 A study on the effect of sample geometry on spatial resolution of the same ATR-FTIR imaging system used in this study has been conducted by Everall et al. 25 for convex solid objects. They determined that this imaging system with a calculated NA of $2 was underestimating the size of 20-140 mm objects and approximating solid spheres, of varying dimensions, to be the same size ($30 to 35 mm). The authors postulated that this was due to the shallow evanescent wave penetration (eqn (1)) and blurring caused by the nite spatial resolution. Our experiment however is different when compared to such a case, as the sample is a rather so solid which, with a gentle anvil pressure, provides a at central area that is quite large ($100 mm) and the convex shape only occurs at the edges. As the main purpose of this paper is not to estimate the real size of the microparticle, but to compare the output of univariate, hard and so multivariate tools, we are more interested in seeing how the different data analysis approaches impact upon measured spatial resolution and the sharpness of the interfaces. The measured spatial resolution for each of the data analysis methods is summarised in Table 1 and the spectra along the grey arrow on all spectra are shown in Fig. S4. † Table 1 shows that the NLCF approach improves the measured spatial resolution of a mid-IR image when there is a large discrepancy between the refractive index in adjoining pixels, i.e. an air/polymer interface. Lasch and Naumann have shown the application of Fourier self-deconvolution to IR microspectroscopy images has led to an improvement in the spatial resolution by improving the spectral resolution at each pixel. 45 The NLCF approach we have used here, is also acting in a spectral resolution enhancement manner and the ability, when used in a supervised manner, to discriminate between species in adjoining pixels more readily than the peak height, MCR-ALS and Gaussian peak tting approaches that limits the blurring effects at interfaces.
The comparison between the images generated using NLCF and the images generated using traditional Gaussian tting is an interesting one. The results are showing that NLCF facilitated an improvement in spatial resolution compared to a Gaussian t. This is likely to be mainly due to the optimisation algorithm (trust region) falling into local minima as a result of a lack of change between consecutive iterations due to the limited number of parameters in the Gaussian tting protocol not allowing an improvement in t, due to variations in the symmetry of bands between pixels. In order to make the algorithm robust for not only this but different data sets and as we know that IR bands are asymmetric, we chose to use the Pearson IV function despite the delay in computation time. The merits of using this function compared to Gaussian and or Lorentzian functions is also discussed in the text and in ref. 38 and 39.  When the refractive index change at an interface is small, such as when a PLGA microparticle is surrounded by water, the blurring of the interface is reduced and the spatial resolution determined by this step-edge approach is comparable between the each of the analysis approaches. It should be noted that for a sharp interface, such as the USAF (1951 1X 38257) target sample described in ref. 44, we have calculated the spatial resolution of the system to be $18 mm using the step edge method.

Image comparison
As indicated previously FTIR imaging is a powerful tool that facilitates the collection of spatially resolved chemically relevant data from samples in real time. By conducting such measurements in situ, we can minimise the perturbations to the system and follow chemical changes from the same region of the same sample. One of the challenges presented is the efficient analysis of the huge data matrices that are generated during such experiments (4096 spectra at each time point). To explore the relative merits of a number of different analysis approaches (univariate peak height image plotting, so and hard multivariate modelling), we have taken single PLGA75/25 microparticle and followed its interaction with water as a function of time at 70 C. The experiment is setup in such a way that the interaction between the particle and water will only occur at the interfaces that we are monitoring. Fig. 4(a) shows ve false colour images obtained using the univariate peak height method and a spectrum taken from close to the 'dry' polymer/hydrated polymer interface in the t ¼ 0 image. From this spectrum we are able to see both of the bands used to determine the polymer distribution (the ester carbonyl at $1745 cm À1 ) and the water distribution (the OH bending mode $1635 cm À1 ). The two images at the extreme le show the distribution of polymer (top) and water (bottom), within the ATR eld of view, immediately aer the experiment was started. Even at this short time (data collection was $5 minutes) there is evidence of an interface layer of hydrated PLGA around the particle, with an apparent concentration gradient from the particle centre outwards towards the aqueous media. As the contact time with water increases, a number of phenomena occur. Firstly the microparticle (dened by the red zone) initially appears to increase in size (2 h), which is indicative of swelling and the hydrated PLGA layer (yellow) becomes thicker. Images collected at times exceeding 2 hours show the particle decreasing in size and the boundaries of the particle becoming less well dened. The complementary images pertaining to water concentration show an inverse relationship, as one would expect in a binary system. The changes to the interface layer as a function of time will be discussed later. It is clear however, that the overlap between the bands used to determine the distribution of the two components in this system must inuence the sharpness of this image and increase the magnitude of the measured interface. Fig. 4(b) shows the equivalent false colour infrared images for the same data described in Fig. 4(a), this time, obtained using a MCR-ALS approach and 2 factors. The spectral features shown to the le hand side of this gure are the 'pure component spectra' for PLGA and water generated using this so modelling method and are oen referred to as factors. This data is in general agreement with the ndings from the peak height measurement approach shown in Fig. 4(a). Closer inspection indicates that the interfaces in this set of images are blurred as were those observed in the univariate data set (Fig. 4(a)) and less sharp than those obtained using the hard modelling approach (Fig. 4(c)). But this data set does exhibit improved S/N compared to the univariate data ( Fig. 4(a)). The blurring of the interfaces is the result of the 'pure component spectrum' representing PLGA still displaying a feature at $1635 cm À1 associated with water and due to the fact that the MCR factor for PLGA has xed peak centres and band widths meaning they approximate rather than exactly replicate the spectrum at each pixel. The improvement in the S/N is the result of the elimination of water vapour in the pure spectral factors combined with the fact that the signal comes from many more spectral data points compared to the peak height data. Fig. 4(c) shows 5 false colour images equivalent to those described in Fig. 4(a) and (b). However, this time they were obtained using the nonlinear curve-tting approach and are the result of the summation of the peaks, generated during the tting process, that have been assigned to PLGA (upper row) and water (lower row). The data to the le of these images shows the 11 component peaks used to t the spectrum, from the same pixel used to obtain the peak height and MCR-ALS data within the PLGA/water interface. The dotted line denotes the synthetic spectrum generated from the combination of the tted peaks, which matches the real spectrum at that pixel.
The upper set of images shows the distribution of PLGA determined using all of the tted peaks except the peak with a maximum at 1635 cm À1 which is used to obtain the distribution of the water within the ATR eld of view. These images are, in general, in agreement with the data shown in Fig. 4(a) in that they indicate the formation of a hydrated region around a dry PLGA particle that increases in thickness over the rst 2 hours and that this occurs concurrently with particle swelling. Closer inspection and comparison with the data in Fig. 4(b) indicates that the interfaces and boundaries in this set of images are much sharper and the data exhibits less noise, i.e. there is less variation in colour intensity between equivalent pixels. This (apparent) improvement in resolution and S/N in each image is achieved by the elimination of contributions from overlapping features at each pixel such as other chemical species, instrument noise and atmospheric water vapour. Another contributing factor to the broadening of the interfaces in the MCR-ALS images when compared with those generated using NLCF, is the fact that the NLCF peak centres are optimised for each peak within each pixel, whereas the MCR-ALS images used a xed factor, with xed peak centres and band widths. It is likely that during ALS optimisation, then a linear combination of the water and PLGA factors may give a better mathematical t in some of the interface pixels resulting in a less well dened image, whereas the NLCF approach is better able to discriminate between water and PLGA. This does come at a considerable time and convenience penalty. The images shown in Fig. 4(a) can be obtained in seconds, whilst the data shown in Fig. 4(c) takes approximately 5 hours to generate using a PC with an Intel® Core™ i7-2620M CPU @ 2.7 GHz and 8 GB of RAM. In many applications this approach may not be feasible due to a number of considerations such as time, CPU availability and, more importantly, spectral data which is too challenging to t due to a lack of knowledge of the species within that system Fortunately, the authors have amassed an understanding of the components within this system, i.e. water and PLGA thus generating condence in the peak assignments. The MCR-ALS approach facilitates the collection of false colour images in a few minutes and thus offers an attractive/acceptable compromise between the slow but accurate hard modelling methodology and the rapid univariate approaches.

Interface analysis
The generation of false colour images, from the mid-infrared imaging dataset highlighted here, facilitates the rapid assimilation of trends in physical processes such as particle swelling, particle shrinkage, hydration layer formation etc. but cannot readily be used to obtain quantitative information about such processes. This can be problematic when a particle is not uniform in shape and oen the dimensions are estimated by assuming a particular geometry (circle, square etc.) that may not be appropriate. To compare the quality of the output generated using the NLCF approach with standard image generation strategies (peak heights and MCR-ALS) we have compared the ndings along the centre line across the particle as a function of time. To facilitate comparison between the data analysis types we have dened two parameters; the full width at half maximum height of the normalised particle prole (A) and the dimensions of the right hand side interface (B). Their derivation is shown in the ESI, Fig. S5. † Fig. 5(a) and (b) show quantitative data extracted from across the centre line of the generated images shown in Fig. 4(a). Fig. 5(a) shows the evolution of the intensity of the water peak at each pixel across the image as a function of time and Fig. 5(b) shows the associated normalised plot of the full width at half height of the parameter 'A' for the images generated using peak heights.
From Fig. 5(a) it is possible to observe that the overall water concentration across the particle increases as a function of time as one might reasonably expect; rising from an intensity $15% of its maximum value at the local minimum at t ¼ 0 h to a value of $75% of its maximum value at the local minimum at t ¼ 9 h. The prole at t ¼ 0, which is $5 minutes aer the particle has been subjected to water, indicates that the initial ingress of water into the particle is rapid, most likely due to the porous nature of the scCO 2 processed starting material. The shape of the prole initially appears to be fairly uniform and becomes less so as time progresses which may reect the irregular shape of the particles. It is likely that initial compression onto the ATR crystal may make the analysed surface uniform (the ATR experiment collects data from the rst 2-10 mm in direct contact with the crystal) but as the particle hydrates, swells and hydrolyses, the signal obtained via the ATR crystal will depend on the volume of the particle directly above the evanescent eld and how it swells and or moves, potentially resulting in a loss of uniformity. Fig. 5(b) shows the complementary data to that in Fig. 5(a) relating to the intensity of the polymer particle extracted from across the centre line of the peak height generated images. As the particle swells, the concentration of polymer measured within any given pixel will decrease and the concentration of water within that same pixel will increase. Therefore the intensity of the polymer peak (which should be the inverse of the water peak in this binary system) would provide an indication of the degree of swelling in the z-direction. Instead, we are using the change in width of the normalised intensity of the peak height of polymer peak, as a function of time, to provide an indication of the degree of swelling in the y-direction. We have chosen to normalise this data as it facilitates the observation of the change in full width at half height maximum (FWHHM) better than the equivalent data with the non-normalised y values, which decrease over time. For the peak height derived images we clearly see the width of this peak increasing as a function of time and this is plotted in Fig. 5(c).
When water is introduced into our system, there exist domains where we only measure water, others where polymer is the dominant signal and others where we observe a clear mixture of water and polymer; a hydrated zone. Determining the exact point where each domain ends and another domain begins is somewhat arbitrary, but some form of denition is necessary if we are to quantitatively compare data extracted from images generated using different approaches. As described above we have dened a hydrated zone 'B' where both the water intensity and the polymer band intensity are below a certain threshold (10% of the maximum value). Fig. 5(d) shows the plot of the B zone for the peak height derived images as a function of time. The size of this zone increases quite dramatically over the course of this experiment, with dimensions around 70 mm at t ¼ 0 and expanding to 140 mm at t ¼ 9 h. This increase in thickness of the outer hydration layer is an interesting nding and in broad agreement with confocal uorescence images generated by Bajwa et al. 46 of hydrating HPMC tablets. In the HPMC system an outer hydration layer of $100 mm increasing to 200 mm was measured over the course of a 'wetting' experiment. Clearly the timescales are different between the two systems due to the inherently different hydrophilicities, but nonetheless this adds credence to the nature of our measurements and our denition of B. Fig. 6(a)-(d) show data comparable to that presented in Fig. 5, this time derived from the MCR-ALS generated images shown in Fig. 4(b). Fig. 6(a) and (b) show the evolution of the water intensity and the normalised plot of the polymer factor intensity respectively. The shape and the intensities of the proles in these gures are similar to those derived from the peak height measurements, but there was less noise in the MCR derived data (most evident in Fig. 6(a)), particularly at longer time points, where the band intensities of the polymer peaks are quite low. This improvement of S/N occurs because the peak height data is derived from a single point and the MCR-ALS data is derived from a large number of data points. There is perhaps some evidence of the polymer band intensity prole being less sharp at longer time points and at its maximum, which could be related to the contribution of water within the extracted pure factor associated with the polymer (MCR factor 2; Fig. 4(b)) and will also be a function of the xed lineshape of factor with its associated peak maxima and minima which will not be able to exactly match the spectrum at each pixel.
Both the B values (Fig. 6(c)) and the full width at half maximum height values (Fig. 6(d)) as a function of time are very similar to those shown in Fig. 5(c) and (d). This indicates that there is perhaps no signicant improvement in the quality of output obtained for this system when performing an MCR-ALS analysis on the data when compared to the more rapid peak height approach. Of course the MCR-ALS method can be used without any prior knowledge of the system; therefore there is no need to identify a peak specic to each component within it, which could be advantageous in some instances. Fig. 7(a)-(d) show data comparable to that shown in Fig. 5  and 6, this time derived from the NLCF generated images (Fig. 4(c)). Fig. 7(a) shows the intensity of the curve tted water band as a function of time. Both the intensities and width of these proles are somewhat different to those observed in those derived from the peak height ( Fig. 5(a)) and MCR-ALS ( Fig. 6(a)). Firstly the intensities are generally higher than those observed for the data derived using the other two approaches, this is more pronounced at short times, with the values at t ¼ 0 being approximately 30% of their nal intensity (cf. 15% for both peak heights and MCR-ALS). Some explanation is found in the consideration of the factors/bands used to generate the initial images from which these line proles were generated. In the case of the peak height data, it is clear that the vector normalisation and baseline correction has enhanced S/N of the images (compare Fig. 4(a) with the data in ESI, Fig. S6 †), but it is entirely feasible that this will exert some inuence on the intensity values generated for each spectrum within a given pixel.
In the case of the MCR-ALS the extracted pure factor for the polymer contains a contribution from water ( Fig. 4(b)) and therefore when the scores at each pixel for the pure water factor are calculated then they will be underestimated. As the NLCF approach is able to generate both a pure water signal and a pure polymer signal free from interference, we anticipate that the intensities presented in these proles will be more likely to match the true concentration prole. Fig. 7(b) is also somewhat different to the analogous peak height and MCR-ALS data; in that it is narrower and different in shape. Once more it is the convolution of the water and polymer bands in both the peak height and MCR-ALS spectra that contributes to the broadening of this prole, relative to the NLCF data. Fig. 7(c) shows the increase in the B zone dimension as a function of time. The values plotted here are lower than those obtained via both the peak height and MCR-ALS approaches and the error determined at each time point is somewhat lower than the corresponding values for the other two methods. This is a clear indication of the ability of the NLCF procedures to eliminate noise from the processed images. Fig. 7(d) shows the FWHM of the normalised polymer band prole generated from the NLCF images ( Fig. 4(c)) and this also shows a reduction in width in comparison with the analogous data generated using the other two methods. The reduction in size of the B zone and the FWHM in comparison with the data generated using peak heights and MCR-ALS is a reection of the ability of the NLCF method to discriminate between the water and polymer contributions at each pixel, reducing the blurring effect of convoluted spectra.

Degradation rate calculation
Coupling the chemical selectivity of infrared spectroscopy with the regional selectivity of an FPA detector facilitates remarkable insight into a wide range of processes. PLGA microparticles can be used as sustained delivery vehicles, where the rate of hydrolysis will control the release rate. Parameters that govern the hydrolytic degradation of PLGA include molecular weight, structure and morphology and PLGA degradation dynamics. FTIR spectroscopy has routinely been used to follow hydrolysis kinetics, but we believe this is the rst time that measurements have been undertaken on single microparticles in this manner. Two infrared bands have been observed at $1452 cm À1 and $1424 cm À1 that correspond to the antisymmetric bending of CH 3 from the lactic acid units and the symmetric bending of CH 2 from the glycolic acid units of the PLGA polymer. The relative intensities of these two bands can be used to estimate the relative quantity of glycolic and lactic acid units present in the polymer and has been used to determine the rate of hydrolysis of the two co-polymer segments within the same experiment. Work by Vey et al. 47 has shown that the lactic acid units hydrolyse $1.3 times slower than the glycolic acid units. Unlike large-sized (a few mm) PLA/GA polymer devices, microspheres less than 300 mm in diameter have been shown to undergo homogeneous degradation with the rate of degradation of the core being equivalent to that at the surface. 48,49 Therefore rate constants from different regions within a microparticle would be expected to give the same calculated rate. Fig. 8 shows typical single pixel spectra ( Fig. 8(a)), spectra resulting from the binning of 5 Â 5 pixels ( Fig. 8(b)), the resultant peak ts from a single pixel (Fig. 8(c)) and the result of peak ts from the binning of 5 Â 5 pixels (Fig. 8(d)) which have been used to calculate the rate constants, shown in ESI, Fig. S6. † It should be noted that due to a deviation from a linear relationship between the ln(intensity) versus time plot, no k values for the MCR data could be determined.
To compare the relative merit of each of the data analysis approaches used in this study, we have calculated the hydrolysis rates for both the glycolic and lactic blocks independently within the same experiment. MCR-ALS was unable to provide pure component spectra for both the lactic and glycolic segments of PLGA, probably due to the nature of the iterative extraction process. MCR-ALS relies on variance within spectral data sets to extract pure component factors and as the ratio between the glycolic and lactic units during hydrolysis is constant throughout the experiment (i.e. the data is co-linear) the algorithm does not detect any variance. Consequently a single 'PLGA' pure component is generated. Hydrolysis rate constants for the lactic units (k L ) and the glycolic units (k G ) were calculated using both the peak height data and the NLCF data from the logarithmic plot of peak intensity versus degradation time (see ESI data S7 †) using the band $1452 cm À1 for the lactic units and that at $1424 cm À1 for the glycolic units, using the spectra shown in Fig. 8. First order kinetics as dened in eqn (10) and (11) were assumed; for the lactic unit and, for the glycolic unit.
Rate constants were calculated by selecting the spectrum (or extracted factor) at 3 random pixels within each image. For comparison, the same positions were used for each of the peak height, MCR-ALS and NLCF image analysis approaches and the mean of the three values are shown in Table 2.
Rate constants were also determined by binning the spectra/ factor score from 3 areas of 5 Â 5 pixels from random regions within each image. Once more the same regions were used for each of the peak height, MCR-ALS and NLCF image analysis approaches and the calculated rate data (where possible) are shown in Table 2 and are the result of the mean of three values.
Tracy et al. 48 studied the degradation of poly(lactide-co-glycolide) microspheres in vivo and in vitro and determined degradation rate constants by measuring the polymer molecular weight as a function of time by gel-permeation chromatography. They found that the in vivo degradation rate was higher than in vitro degradation one. Their calculations of in vivo rate constant by GPC analysis for ester capped and uncapped (-COOH) PLGA50:50 microspheres were 0.033 AE 0.006 per day and 0.13 AE 0.05 respectively. Considering that the PLGA microparticle studied here was scCO 2 processed therefore very porous, and a calibration (i.e. %GA or %LA versus IR absorbance) was not considered as in ref. 47, due to having one polymer composition, the k values calculated are in reasonable agreement with each other and with the values determined by Tracy et al. 48 and, as one might anticipate, the error obtained when calculating rate constants by binning a number of spectra Fig. 8 Typical single pixel and averaged (over 5 Â 5 pixels) spectra of pre-treated raw ((a) and (b)) and non-linear curve fitted ((c) and (d)) data. Spectra from bottom to top were obtained at t ¼ 0, t ¼ 2 h, t ¼ 5 h, t ¼ 7 h and t ¼ 9 h, respectively. Table 2 List of degradation rate constants (day :1 ) calculated for triplicates of averaged 5 Â 5 pixels and single pixel spectra using peak height values and curve fitted area values over the course of the 9 h hydrolysis experiment. Errors quoted are the standard deviation of 3 measurements taken from 3 different regions of the same particle  is somewhat lower than that obtained from a single pixel. The errors are larger for the peak height derived calculations than those from the NLCF measurements. It is also clear that the ratio of the rate constants calculated for the lactic and glycolic groups is also comparable with that determined by Vey et al.; 47 we have determined the ratio of rate constants to be 1.2 (peak height binned pixels), 1.7 (peak height single pixels), 1.4 (NLCF binned pixels) and 1.3 (NLCF single pixels). Interestingly the error determined for the k calculations for the MCR-ALS processed images seemed to be independent of the number of pixels used to determine them. It is unclear if this nding is real or an anomaly of the pixels chosen to make the measurements.

Conclusions
A new algorithm for nonlinear curve tting (NLCF) has been developed and applied to the analysis of mid-IR images of a single PLGA microparticle undergoing hydrolysis at 70 C for the rst time. The supervised NLCF approach has been shown to have several advantages over traditional peak height measurements and a commonly applied multivariate tool; MCR-ALS. Firstly the application of NLCF routines to such data has been shown to enhance the spatial resolution within a sample with (a) overlapping spectral signals and (b) containing interfaces with large discrepancies in the refractive index (i.e. air/polymer). Secondly it has been shown to improve S/N and sharpen features such as interfaces in processed images, due to its ability to discriminate between different species in a mixture, this is particularly pronounced when one compares this approach to standard peak height measurements. Thirdly as the approach does not appear to be greatly inuenced by colinearity, unlike MCR-ALS, supervised NLCF can be used to extract chemical information from species changing at the same ratio during a kinetic process such as hydrolysis. Finally the high S/N at each pixel readily facilitates the calculation of rate constants from a single pixel with a low error when compared to traditional peak height approaches. All these advantages come at a signicant time penalty; the NLCF algorithm described takes $5 hours to extract information from a single mid-IR image containing 4096 spectra, this compares to $2 seconds for peak height analysis and $1 minute for MCR-ALS on the same image using the same PC.
With the help of fast developing computing hardware power, the supervised NLCF method developed here will be a useful IR imaging analysis tool providing high resolution images and quantitative analysis for many more cases particularly where hard modelling is the only option such as deconvoluting protein spectra when searching for changes in secondary structure.