Article pubs.acs.org/crystal
Image Analytical Approach for Needle-Shaped Crystal Counting and Length Estimation Jian X. Wu,†,# Sergey V. Kucheryavskiy,‡,# Linda G. Jensen,† Thomas Rades,† Anette Müllertz,† and Jukka Rantanen*,† †
Department of Pharmacy, Faculty of Health and Medical Sciences, University of Copenhagen, Universitetsparken 2, 2100 Copenhagen, Denmark ‡ Section of Chemical Engineering, Department of Chemistry and Bioscience, The Faculty of Engineering and Science, Aalborg University, 9100 Aalborg, Denmark ABSTRACT: Estimation of nucleation and crystal growth rates from microscopic information is of critical importance. This can be an especially challenging task if needle growth of crystals is observed. To address this challenge, an image analytical method for counting of needle-shaped crystals and estimating their length is presented. Since the current algorithm has a number of parameters that need to be optimized, a combination of simulation of needle crystal growth and Design of Experiments was applied to identify the optimal parameter settings for the algorithm. The algorithm was validated for its accuracy in different scenarios of simulated needle crystallization, and subsequently, the algorithm was applied to study the influence of additive on antisolvent crystallization. The developed algorithm is robust for quantifying heavily intersecting needle crystals in optical microscopy images, and has the potential for being applied for real-time crystallization monitoring and control.
1. INTRODUCTION Measurement of crystal morphology is of critical importance, for example, for the fine chemicals and pharmaceutical industries. Particle size and shape distributions can affect downstream processability, as well as other performance characteristics such as the dissolution rate of pharmaceutical products. Several fast and nondestructive analytical techniques1 for monitoring2−5 and control of crystallization process have been reported, including spectroscopy,6,7 laser diffraction,8,9 and imaging10−13 techniques. When these techniques are used in a complementary fashion, insight into several phenomena, such as onset of crystallization, crystal growth rate, solution phase concentration, and level of supersaturation, as well as solid form composition can be investigated.11,14 Although the prospect of achieving crystallization control is promising, the critical challenge in the pharmaceutical sector is the limited amount of sample available in the early phase of development,15 and therefore, small-scale vision based screening platforms have been proposed to study factors affecting crystallization.16−18 This would be of special interest not only for process development, but also for the early drug development phase where different additives capable of delaying the crystallization are utilized as a part of pharmaceutical product development.19,20 Several studies have demonstrated the applicability of nondestructive real-time monitoring of nucleation related phenomena and crystal growth using microscopy based © XXXX American Chemical Society
monitoring technique together with image analytical routines.12,21−23 When compared to subjective assessment of the microscopy images,16,24,25 the image analytical techniques are objective allowing fast extraction of numerical information from the images and hence these methods are less dependent on the assessment from an expert panel.14 Although promising, earlier studies11,22,23 have revealed that the precision of the developed image analytical method is dependent on the crystal density as well as on the applied image analytical routine. As crystal density within the field of view increases, and eventually crystals are touching each other, multiple crystals will ultimately appear as a single object in the image leading to poor estimation of the amount of crystals as well as compromised size estimation.11,22 This problem is more pronounced for needle-shaped crystals, where the tendency of needles intersecting each other will affect the conventional image counting algorithms,22,26 making them inapplicable for counting the exact number of needles in the image. Although effort has been spent solving this problem as in the study by Larsen et al.,23 the proposed method is more complex from an image analytical point view as well as computationally more expensive. Hence there is a need to develop an image analytical routine that is simpler to implement, and that allows counting of needle-shaped crystals Received: May 26, 2015 Revised: August 30, 2015
A
DOI: 10.1021/acs.cgd.5b00720 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
with slope a and intercept b in a Cartesian coordinate system with origin in the center of an image plane. Although eq 2
as well as quantifying their length to estimate the extent of nucleation and crystal growth. In order to extract information related to line features from the image, line detection algorithms such as Radon transformation26 are well-established image analytical methods. Although Radon transformation has been applied in other scientific areas such as estimating X-ray projections and tomographic reconstruction,27 its applicability as well as robustness for counting and estimation of needle-shaped crystals remains to be explored. In an earlier study, it was demonstrated that computer simulation of reference images is a convenient approach to estimate the accuracy of the output from an image analytical routine,22 and furthermore, the importance of critical examination of the validity of the developed image analytical routine prior to its application has been stressed.14 Since many mathematical parameters exist for the Radon transformation routine, there is a need to estimate the combination of these parameters, so that the developed routine has the optimal accuracy in relation to its application for estimating the nucleation and crystal growth rate for needleshaped crystals. The aim of present study was to create, optimize, and validate a Radon transformation based image analytical method for real-time monitoring of needle-shaped crystals (count and length). To demonstrate the applicability of the developed method, the influence of a polymeric additive polyvinylpyrrolidone (PVP) on the carbamazepine antisolvent crystallization has been exemplified.
yi = axi + b
can be used to describe this line, in practice when the line approaches vertical, the slope value a becomes infinitely large. Hence from a computational perspective, it is practical to express a line using the mathematical expression in eq 3,26 which can be derived from eq 2
x cos(θ) + y sin(θ) = ρ
Figure 1. Illustration of a black line expressed using the angle θ and ρ. The Radon transformation in its general form is defined on a space of straight lines L in R2 as a function of two argumentsslope of the lines (θ) and their right angle distance from origin to the line (ρ). The transformation maps any differentiable function f(x,y) from the R2 space to the function f ′(θ,ρ) in the space L according to eq 426
2.1. Materials. Hydroclorid acid (37%) was purchased from VWR (Søborg, Denmark) and carbamazepine was purchased from Fagron (Copenhagen, Denmark). PVP was obtained from BASF (Kollidon 90F, Ludwigshafen, Germany). All of the above materials were of Ph. Eur. grade. N,N-Dimethylacetamide was of analytical grade and purchased from Sigma-Aldrich (Brøndby, Denmark). 2.2. Carbamazepine Antisolvent Crystallization. Antisolvent crystallization without the influence of PVP was performed by dissolving carbamazepine in N,N-dimethylacetamide and 0.5 μL of the solution was placed on a microscope slide. The crystallization was subsequently performed by adding 19.5 μL 0.1 M hydrochloric acid diluted with aqua to the carbamazepine solution, and placing a cover glass (22 × 22 mm2) on top. In order to investigate the influence of polymer on the carbamazepine crystallization, another set of experiments was performed, where 25% (w/w) PVP was dissolved together with carbamazepine in N,N-dimethylacetamide. The crystallization was monitored using a cross-polarized light microscope (Axiolab, Carl Zeiss, Göttingen, Germany), and a video was captured at 2.67 frames per second for 15 min using a camera (Moticam 10MP, Motic, Wetzlar, Germany) interfaced with the software Motic Images (Motic Images Plus 2.0, Motic, Wetzlar, Germany). Each frame has dimension 1832 μm × 1374 μm. The objective used in the microscope had a 5× magnification and 0.12 numerical aperture. All experiments were performed as triplicates. To compare the kinetics of the observed phenomena, the needle count and average needle length profiles were fitted with a logarithmic growth model using a Levenberg−Marquardt algorithm.28 The fitted model is presented in eq 1:22 α 1 + exp(− β × (t − t50 ))
(3)
where the angle θ and parameter ρ are related to the slope and the right angle distance from origin to the line (Figure 1), respectively.
2. MATERIALS AND METHODS
X=
(2)
+∞
f ′(θ , ρ) =
+∞
∫−∞ ∫−∞
f (x , y)δ(x cos(θ) + y sin(θ) − ρ) dx dy (4)
This can be specified for the case where the arguments x, y as well as θ and ρ are discrete, which makes it applicable for image analytical application, since pixel values in the images are discrete elements. To demonstrate the principle behind Radon transformation in line detection, first an image plane is defined as a Cartesian coordinate system (x, y), with origin in the center of the plane as previously described. Consider a line object in the image with the function f(x,y), where f represents intensity of a pixel with coordinates (x, y). Not to be confused with the line of interest to be detected in the image just mentioned, consider a set of n lines, crossing the origin of the coordinate system each forming an angle θ with the abscissa, given as θ = {θ1, θ2, θ3, ..., θn}, where θi is an angle between the ith line and the abscissa. On every line a set of equidistant points r = {−rm, −rm−1, ..., −r2, −r1, 0, r1, r2, ..., rm−1, rm} is selected, where rj is a distance to the jth point from the origin on the line. The discrete Radon transformation f ′(θ,r) of the object f(x,y) for arguments (θi,rj) can be calculated as a sum of intensities of pixels from f, whose projections to the line, defined by θi, lie within (rj − d/2, rj + d/2). Here d = |rj − rj−1| is a resolution parameter for r. The projection of a pixel with coordinates x and y to the line is according to eq 3. Therefore, all pixels, whose projections meet the condition in eq 5 |x cos(θi) + y sin(θi) − rj| < d /2
(5)
will be counted. Figure 2 illustrates a schematic example of the working principle behind the Radon transformation for two segments with different slopesone orthogonal to the projection line and crossing it at the rj, defined by θi and one is not. Each segment consists of four pixels represented as blue points in the Figure 2. Assuming all pixels on the image object have an intensity equal to one, then f ′(θi,rj) = 3. The inset in Figure 2 presents the result of the transformation for all points r where there is at least one projected pixel. By repeating this
(1)
where X, α, β, t, t50 are the fitted response, the response value at equilibrium, growth rate constant, time, and time for 50% converted, respectively. 2.3. Radon Transformation. To set the scene for introducing the Radon transformation in detecting line features in an image, first suppose a line is formed from N discrete coordinates (x1,y1)...(xN,yN) B
DOI: 10.1021/acs.cgd.5b00720 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
required. The algorithm presented below is written in an in-house developed Matlab (v 8.2 MathWorks, Nattick) script together with the Image Processing Toolbox (v 8.3, MathWorks, Nattick) and Statistics Toolbox (v 9.1, MathWorks, Nattick). Table 1 presents the main steps of the proposed algorithm. When the raw image data is a grayscale image, the aim of the first part of the algorithm (steps 1−5 in Table 1) consists of creating a binary image, where the noise, e.g., from the detectors and artifacts in the image, is minimized. After having preprocessed according to steps 1−5 in Table 1, the remaining features in the binary image are mainly related to the thin intersecting needles. Since some of the preprocessing steps consist of mathematical morphological operations generally applied in image analysis, they will be discussed briefly in the present work, and details can be found in the work by Gonzalez and Woods.26 The algorithm steps 7−9 include the actual Radon transformation previously described, with the aim of obtaining the image mapping in (θ, r) space. The final step 10 (Table 1) consists of reconstructing the line, based on both the identified peaks (θ, r) pairs and the binary pixels from the output binary image in step 5 (Table 1). This step is performed by identifying all pixels from the binary image that satisfy the (θ, r) criteria. Using these pixels as starting point, other pixels are assigned as belonging to the same line based on distance criteria. It has been found that the computational speed can be increased by splitting the main image into several subimages (step 6), and process each of them independently through steps 7−10. In the following, each of the steps in Table 1 will be discussed together with illustrated examples demonstrating the effect of the preprocessing. In step 1, the original grayscale image (Figure 4a) is binarized by applying a hard threshold, so that all pixel values below and above the threshold are set to 0 and 1, respectively. In the obtained binary image (Figure 4b), it can be observed that the needles in the foreground, which are in focus, have been segmented, while the needles in the background with lower intensity values have been left out from further analysis. Followed by morphological opening and closing operations,14,26 it is observed that many small artifacts have been removed from the image (the green boxes in images belonging to steps 1−2 in Figure 4b and Figure 4c), which improves the accuracy of the subsequent Radon transformation by decreasing the number of pixels that do not lie on a line. In order to obtain the main linear features, the lines are thinned using a mathematical operation known as
Figure 2. Scheme and result of Radon transformation of an image object f for slope θi. The illustration demonstrates six pixels (illustrated with blue dots) forming two orthogonal lines (illustrated with black lines labeled as 1 and 2) being projected onto one of the n lines crossing the origin with slope θi. The inset illustration is an example of the projection results, where both lines 1 and 2 are projected. Because line 1 is orthogonal to the line crossing the origin, its projection will results in a local maxima at rj, while this will not be the case for line 2. procedure for all angle values a full representation of the function in Radon space can be achieved. Thereby, any digital image can be mapped from the Cartesian space (x,y) to the Radon space (θ,r). Hence intersecting line features in an image can easily be identified as separate lines when mapped to the Radon space followed by estimation of the (θ,r) as demonstrated in Figure 3. This basic principle forms the basis for building a needle-shaped crystal detection algorithm described in further detail below. 2.4. Needle-Shaped Crystal Detection Algorithm. Although the principle behind Radon transformation as described previously appears simple at first glance, in order for this algorithm to be reliable and accurate in detecting needle-shaped crystals, a number of image preprocessing steps and optimization of the algorithm’s parameters are
Figure 3. Radon transformation in separating intersecting lines. Left: Three intersecting lines that if counted using ordinary counting algorithm will result in one object. Right: When the original image (left) was mapped to the Radon space, each intersecting lines result in a local maximum which can be counted. The lines and their corresponding local maxima in the Radon space are labeled. Origin of the left image is located at the image center. C
DOI: 10.1021/acs.cgd.5b00720 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
Table 1. Overview of the Image Analytical Steps and the Algorithm Parameters for the Needle Detection Algorithm step
inputs
1 2
grayscale image output step 1
thresholding morphological opening followed by closing
operation
tunning parameters
3 4
output step 2 output step 3
skeletonization pruning
5
output step 4
dilation
6
output step 5
7 8
output step 6 output step 7
9 10
output step 8 output step 9
split the image into sub images, and process each of them independently radon transformation Gaussian filtering followed by subtraction of output in step 6 with filtered matrix thresholding, peak identification in Radon space line reconstruction in original image
threshold value disk shaped morphological structuring element -
skeletonization26 (step 3 in Figure 4d). Although useful in extracting the main linear features, the skeletonization process can introduce artifacts in the image seen as small branches on the main needle (highlighted in the green box in the image belonging to step 3 in Figure 4d). These branches can be minimized (highlighted in the green box in the image belonging to step 4 in Figure 4e) using a pruning step.26 Although the Radon transformation can be performed on the pruned image (Figure 4e), in practice, it is found that the prediction accuracy is increased when the linear features on the pruned image have been thickened using a dilation step26 (step 5 in Figure 4f). This is because the number of pixels that are situated along a line is increased by image dilation, hence increasing the signals to noise ratio in the Radon space, making it easier to extract the (θi,ri) local maxima (step 7 in Figure 4g). Prior to extraction of the local (θi,ri) maxima, in order to minimize the high frequency noise in Radon space obtained in step 7, the Radon space is filtered using a rotationally symmetric Gaussian filter of size 80 × 80 pixels and 20 in standard deviation. Hereafter the original Radon space is subtracted from the filtered Radon space, and the local (θi,ri) maxima extracted using a hard threshold. From the identified (θi,ri) pairs the pixels in the binary image lying along the line can be identified by applying eq 4. In some cases, due to the thresholding in the image processing steps, the line might appear to be broken into several segments (observe an example in Figure 4). In order to minimize this problem, the identified pixels are sorted according to their position on the line. This is followed by calculating distances between successive pixels, and if the distance is below a predefined gap value, the line segments will be joined together. Although most of the lines are identified in this step, this algorithm can introduce a number of wrong detected lines, which is due to some pixels although not lying along the line that belongs to a needle crystal, by chance can form false lines, which looks like small “branches” alongside the main line. Such artifacts, e.g., can be caused by the skeletonization step, and can be observed (Figure 4 h). In order to minimize this problem, a minimum line length threshold has been introduced, and detected line lengths below this threshold are removed from the image. This step greatly minimized the problem related to false detection, and improves the accuracy of the algorithm. 2.5. Validation and Optimization of the Needle Detection Algorithm. The above proposed algorithm has several parameters, which may influence its accuracy on the estimated needle counts and lengths. Hence it is important to evaluate the influence of algorithm parameter values on its overall estimation accuracy. In order to estimate the optimal algorithm parameter values, a Design of Experiment (DoE) combined with image simulation approach has been applied. As this approach requires reference values, a set of binary images with predefined number of lines and lengths were simulated and the needle detection algorithm with different parameter value combinations according to the DoE was applied to estimate the needle counts and average length of the reference image data set. Since real-life needle crystals rarely are fully straight; bends were also introduced into these lines using quadratic Bezier’s curves. The
disk-shaped morphological structuring element number of splits θ increment Gaussian filter design threshold value minimum line length, line fill gap
output binary image an image with reduced noise and main objects filled binary skeletonized image with spurs removing spurs from the skeletonized image thickens the needles, so that more binary pixels lie on a line splitted images radon parameter space Radon parameter space with decreased high frequency noise (θ, r) line count and length
shape, orientation, and position of such curves are defined by three control points P1, P2, and P3, located such that they form an isosceles triangle (see Figure 5 for an example). The curve starts at P1 and ends at P3, while P2 defines the curvature extent. The extent of curvature is defined as the ratio between triangle height from the base P1P3 to P2 and half the line segment |P1P3|. The points of the curve are calculated according to eq 6 ⎧ B = (1 − t )2 P + 2(1 − t )tP + t 2P ⎪ x 1x 2x 3x ⎨ ⎪ By = (1 − t )2 P1y + 2(1 − t )tP2y + t 2P3y ⎩
(6)
where t is a continuous parameter between 0 (the calculated point will be at P1) to 1 (where the calculated point will be at P3). Bx and By are the x and y position of the pixels for the bended curve. Figure 5 presents examples of generated quadratic Bezier curves with different location of the control points. The curves were positioned randomly on the simulated images and had a random slope. Images with all combinations of curvature {0%, 3%, 5%}, number of needles {10, 200, 400, 600, 800, 1000, 1500, 2000}, and lengths {20, 50, 100, 150, and 200 μm} were simulated using an in-house-developed Matlab script. Figure 6 presents several examples of simulated images. The parameters from the needle detection algorithm investigated and the levels of variation of the full factorial design are presented in Table 2. The full factorial design consisted of 320 runs and the accuracy of the predicted needle count and mean needle length where calculated according to eq 7
accuracy = 1 −
|ref − estimated| ref
(7)
where “ref” and “estimated” are the reference value of either simulated needle count or average needle length and the estimated value from experiments carried out according to the DoE scheme, respectively.
3. RESULTS AND DISCUSSION 3.1. Optimal Parameters for the Needle Detection Algorithm. The influence of varying the following effects; splitting the image into several subimages, threshold level, minimum line length, and theta step resolution; on the estimated line count and average line length accuracy, is illustrated in the effect plot (Figure 7). For both responses, it was found that splitting the image into several subimages not only increased the data processing speed, but also improved the overall accuracy (Figure 7). The main reason for this is that while, splitting the main image into several subimages, the size of Radon space belonging to the transformed subimages remains the same, the number of peaks in the Radon space is decreased. Hence, it is easier for the peak detection algorithm to identify and extract the (θ,r) maxima. Although the approach D
DOI: 10.1021/acs.cgd.5b00720 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
Figure 4. Example of the needle detection algorithm from carbamazepine crystallization. The preprocessing steps of the image and line detection: (from top to bottom and from left to right) original grayscale image (a), binarization by thresholding (b), opening and closing (c), skeletonization (d), pruning (e), dilation (f), result of Radon transformation (g), detected lines presented as the green lines in (h) using the detected (θ,r) maxima from (g). The green boxes between two successive images highlight the effect of the preprocessing step.
to peak detection in the Radon space from 0.2 to 0.35 increased the accuracy for the needle count estimation. The reason for this observation is that upon an increase in the threshold value, only high intensity peaks mainly belonging to the needles are kept. However, low intensity peaks, which mostly belong to the pixels that by chance are arranged along a line are removed from being detected. While the threshold factor has an influence for count response, variation of this factor within the specified range has no influence on the average needle length response. It was hereafter decided to apply a threshold value of 0.35 for the line count and average line length estimation. For the threshold related to minimum line length detected, it was observed for both responses that an increase in the minimum line length decreased the accuracy (Figure 7). The main reason behind this finding is that initially an increase in the minimum line length mostly removes previously
with subimages was useful for improving the overall accuracy within the evaluated range of numbers of needles and their average length in the simulated images, it should be noted when few needles were present in the image, e.g., at the beginning of crystallization, splitting the image can lead to a situation where part of a needle was present in both subimages, leading to an overestimation and underestimation of needle count and average needle length, respectively. For the subsequent analysis of the experimental needle crystallization, it was decided to split the image into 24 subimages (Figure 7). An increase in the θ resolution from 0.05° to 0.35° did not improve the average length estimation accuracy, and although a slight improvement in the trend for needle count accuracy was observed, this was not significant. It was decided to keep the θ resolution at 0.35° for both responses hereafter. An increase in the threshold value (corresponding to step 9 in Table 1) related E
DOI: 10.1021/acs.cgd.5b00720 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
Table 2. Overview of the Factors, Their Role in the Image Analytical Steps Described in Table 1 and Their Level of Their Variation in DoE factors
step
levels
θ increments Radon parameter space threshold value image split minimum line length fill gap
6 8
0.05, 0.2, 0.35 0.2, 0.35, 0.5 (1, 1), (2, 3), (4, 6) 15, 20, 25 15, 20, 2
9 9
found, it must be emphasized that these optimum parameters are dependent on image size and also the experimental conditions, e.g., the level of objective magnification applied. 3.2. Validation of the Needle Detection Routine. The importance of validating the conditions in which the developed image analytical routine leads to satisfactory results has been highlighted by Wu et al.14,22 Using the optimal parameters for the needle detection routine identified by the present DoE, the accuracy of the estimated average needle length and needle count under the variation of reference needle count and average needle length in the simulated images is presented in Figure 8A and B, respectively. In order to highlight the difference and the importance of understanding the image analytical routine and validate the experimental space before uncritically applying a given image analytical method, a previously reported image analytical method17,22,29,30 has also been tested on the simulated images and its accuracies are presented in Figure 8C and D. As observed in the accuracy map (Figure 8), although both approaches can be used to analyze the images, applying a method, which does not fit for the purpose of
Figure 5. Example of quadratic Bezier curve generated with three control points P1, P2, and P3 forming an isosceles triangle. For the purpose of clear illustration, the extent of curvature is 67% in this case, which means that the ratio between the triangle height from P1P3 to P2 and half line segment of |P1P3| is 0.67.
described false lines created. When this parameter is further increased, both the false created lines as well as true lines belonging to the needles are removed, leading to a reduction in accuracy estimation. Hence for the subsequent analysis of crystallization image data, the optimum minimum line length detected for both responses was set at 15 pixels. Although the optimum parameters for the image analytical routine have been
Figure 6. Example of simulated images with different settings. (a) count 500, average length 100 μm, curvature 0; (b) count 1000, average length 100 μm, curvature 3; (c) count 2000, average length 50 μm, curvature 5; (d) count 3000, average length 200 μm, curvature 3. F
DOI: 10.1021/acs.cgd.5b00720 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
Figure 7. Main effects plot illustrating the influence of variations in the algorithm parameters on the accuracy related to needle count (top) and average needle length estimation (bottom).
Figure 8. Accuracy estimation of the developed needle detection algorithm under the variation of different count and average needle length from simulated images. A and B: Accuracy of the average needle length and count estimation using the needle detection algorithm. C and D: Accuracy of the average needle length and count estimation using the image analytical method previously described.22
is between 60% and 90% (Figure 8A and C). Similarly, a dramatic increase in accuracy when estimating the number of needles in the image has also been found when comparing the previously published method22 with the needle detection routine (Figure 8B and D). This discussion does not invalidate the previously published method. Rather, it emphasizes the importance of understanding the challenges behind an image data set prior to designing the image analytical solution with the
analyzing needle-shaped crystals, leads to poor accuracy estimations for both responses. Comparing the previously reported method22 with the present method in estimating average needle length, it can be seen that when the number of needles in an image exceeds around 200 and the average needle length using the present applied magnification exceeds 50 μm, using the previously reported method leads to almost 0 accuracy, while the accuracy using the needle detection routine G
DOI: 10.1021/acs.cgd.5b00720 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
Figure 9. Representative experimental polarized light micrographs of carbamazepine antisolvent crystallization at the equilibrium. The scale bar corresponds to 100 μm. A and B correspond to carbamazepine antisolvent crystallization without and with polymer (PVP), respectively.
Figure 10. Carbamazepine antisolvent needle crystallization kinetic profiles obtained using the developed needle detection algorithm. Red: with PVP, blue: without PVP. A: Needle count. B: Average needle length.
estimation starts at around 20 μm, which is due to the applied minimum line length in the needle detection algorithm (set at 15 μm)as mentioned previously, in order to minimize the false detection, needles below 15 μm are not counted in the average length estimation. This decision has the drawback that the initial part of the needle growth is not included in the developed image analytical method; however, this problem can practically be solved by increasing the degree of magnification in the polarized light microscope. By following the numerical features extracted from the images over time, kinetic profiles representing nucleation and crystal growth rate were obtained (Figure 10). The fitted α- and β-values for the responses at equilibrium and growth rate constant, respectively, are presented in Table 3. The regression coefficients of fitting were generally above 0.95 indicating good fit to the experimental data. As observed in Table 3, in the presence of PVP in the antisolvent, the needle count (estimate of nucleation extent) is increased at equilibrium as compared with the experiment without PVP in the antisolvent. This observation together with the smaller needles obtained in the
aim of solving a practical problem. Using the optimal parameters for the needle detection algorithm, the computational time to extract count and average needle length of one image frame with 2000 needles with average length of 60 μm is around 5 s on a computer with 2.70 GHz CPU and 8 GB RAM. 3.3. Monitoring of a Model Crystallization System. An example of the outcome of carbamazepine antisolvent crystallization at equilibrium without and with the PVP is presented in Figure 9. It can be observed from the polarized light micrographs that a higher extent of nucleation occurred in the presence of the model polymer in the antisolvent (Figure 9). Furthermore, it can be observed from the polarized light micrograph that when the experiments were carried out with PVP, the needle length in general is smaller as compared with experiments carried out without PVP (Figure 9). The visual impression obtained by examining the polarized light micrograph (Figure 9) was confirmed by examining the count and average needle length profiles obtained by analysis of the images followed by numerical feature extraction (Figure 10). As can be observed in Figure 10B, the average needle length H
DOI: 10.1021/acs.cgd.5b00720 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
Assessment of Recent Process Analytical Technology (PAT) Trends: A Multiauthor Review. Org. Process Res. Dev. 2015, 19, 3−62. (2) Sigdel, M.; Pusey, M. L.; Aygun, R. S. Real-Time Protein Crystallization Image Acquisition and Classification System. Cryst. Growth Des. 2013, 13, 2728−2736. (3) Borchert, C.; Temmel, E.; Eisenschmidt, H.; Lorenz, H.; SeidelMorgenstern, A.; Sundmacher, K. Image-Based in Situ Identification of Face Specific Crystal Growth Rates from Crystal Populations. Cryst. Growth Des. 2014, 14, 952−971. (4) Singh, M. R.; Chakraborty, J.; Nere, N.; Tung, H.-H.; Bordawekar, S.; Ramkrishna, D. Image-Analysis-Based Method for 3D Crystal Morphology Measurement and Polymorph Identification Using Confocal Microscopy. Cryst. Growth Des. 2012, 12, 3735−3748. (5) Roelands, C. P. M.; ter Horst, J. H.; Kramer, H. J. M.; Jansens, P. J. Analysis of Nucleation Rate Measurements in Precipitation Processes. Cryst. Growth Des. 2006, 6, 1380−1392. (6) Heinz, A.; Strachan, C. J.; Gordon, K. C.; Rades, T. Analysis of solid-state transformations of pharmaceutical compounds using vibrational spectroscopy. J. Pharm. Pharmacol. 2009, 61, 971−988. (7) Pataki, H.; Csontos, I.; Nagy, Z. K.; Vajna, B.; Molnar, M.; Katona, L.; Marosi, G. Implementation of Raman Signal Feedback to Perform Controlled Crystallization of Carvedilol. Org. Process Res. Dev. 2013, 17, 493−499. (8) MacLeod, C. S.; Muller, F. L. On the Fracture of Pharmaceutical Needle-Shaped Crystals during Pressure Filtration: Case Studies and Mechanistic Understanding. Org. Process Res. Dev. 2012, 16, 425−434. (9) Guchardi, R.; Frei, M.; John, E.; Kaerger, J. S. Influence of fine lactose and magnesium stearate on low dose dry powder inhaler formulations. Int. J. Pharm. 2008, 348, 10−17. (10) Leyssens, T.; Baudry, C.; Escudero Hernandez, M. L. Optimization of a Crystallization by Online FBRM Analysis of Needle-Shaped Crystals. Org. Process Res. Dev. 2011, 15, 413−426. (11) Caillet, A.; Rivoire, A.; Galvan, J.-M.; Puel, F.; Fevotte, G. Crystallization of Monohydrate Citric Acid. 1. In Situ Monitoring through the Joint Use of Raman Spectroscopy and Image Analysis. Cryst. Growth Des. 2007, 7, 2080−2087. (12) Qu, H.; Louhi-Kultanen, M.; Kallas, J. In-line image analysis on the effects of additives in batch cooling crystallization. J. Cryst. Growth 2006, 289, 286−294. (13) De Anda, J. C.; Wang, X. Z.; Lai, X.; Roberts, K. J.; Jennings, K. H.; Wilkinson, M. J.; Watson, D.; Roberts, D. Real-time product morphology monitoring in crystallization using imaging technique. AIChE J. 2005, 51, 1406−1414. (14) Wu, J.; van den Berg, F.; Rantanen, J.; Rades, T.; Yang, M. Current Advances and Future Trends in Characterizing Poorly Watersoluble Drugs Using Spectroscopic, Imaging and Data Analytical Techniques. Curr. Pharm. Des. 2014, 20, 436−453. (15) Balbach, S.; Korn, C. Pharmaceutical evaluation of early development candidates “the 100 mg-approach. Int. J. Pharm. 2004, 275, 1−12. (16) Van Eerdenbrugh, B.; Taylor, L. S. Small scale screening to determine the ability of different polymers to inhibit drug crystallization upon rapid solvent evaporation. Mol. Pharmaceutics 2010, 7, 1328−37. (17) Xia, D.; Wu, J. X.; Cui, F.; Qu, H.; Rades, T.; Rantanen, J.; Yang, M. Solvent-mediated amorphous-to-crystalline transformation of nitrendipine in amorphous particle suspensions containing polymers. Eur. J. Pharm. Sci. 2012, 46, 446−454. (18) Zhu, Q.; Toth, S. J.; Simpson, G. J.; Hsu, H.-Y.; Taylor, L. S.; Harris, M. T. Crystallization and Dissolution Behavior of Naproxen/ Polyethylene Glycol Solid Dispersions. J. Phys. Chem. B 2013, 117, 1494−1500. (19) Brouwers, J.; Brewster, M. E.; Augustijns, P. Supersaturating drug delivery systems: The answer to solubility-limited oral bioavailability? J. Pharm. Sci. 2009, 98, 2549−2572. (20) Thomas, N.; Holm, R.; Müllertz, A.; Rades, T. In vitro and in vivo performance of novel supersaturated self-nanoemulsifying drug delivery systems (super-SNEDDS). J. Controlled Release 2012, 160, 25−32.
Table 3. Summary of the Fitted Kinetic Parameters Related to Extent of Nucleation and Crystal Growth experiment
count [α]
count [β, s−1]
average length [α, μm]
average length [β, s−1]
No PVP PVP
608 ± 21 1448 ± 446
2.0 ± 0.6 2.0 ± 0.3
55.0 ± 4.6 38.2 ± 3.0
1.4 ± 0.3 0.8 ± 0.2
presence of PVP indicated that the role of PVP is to retard crystal growth by interfering with the fastest growing phase of the needlean observation which is in agreement with previous studies.17,29,31,32 Hence, crystallization of carbamazepine in the presence of PVP is mainly dictated by the nucleation event. Although a difference in the estimate of nucleation rate was not found in the present study (by comparing β-values for needle count in Table 3), the experiments suggested that in the presence of PVP the crystal growth rate was decreased (by comparing β-values for average needle length in Table 3). This observation together with a significant decrease in needle length in the presence of PVP (Table 3) can explain the reason why the number of nuclei is higher in the presence of PVP.
4. CONCLUSION Radon transformation is a useful image analytical method for analyzing the crystallization of needle-shaped crystals. For the algorithm to be practically useful, it is important to preprocess the images in a rational way as well as to optimize the parameters included in the algorithm. The Design of Experiment together with simulation of needles proved to be a robust approach for optimization of the algorithm, and it provided a deeper insight into the performance of the developed algorithm. The performance of the algorithm was demonstrated with antisolvent crystallization of carbamazepine and the effect of additive polymer in the crystallization medium could be investigated.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Phone: +45-35336585. Fax: +45-35336030. Author Contributions #
Jian X. Wu and Sergey V. Kucheryavskiy contributed equally to this work. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS Funding from the Danish Council of Technology and Innovation for the Innovation Consortium NanoMorph (952320/2009) is acknowledged.
■
REFERENCES
(1) Simon, L. L.; Pataki, H.; Marosi, G.; Meemken, F.; Hungerbühler, K.; Baiker, A.; Tummala, S.; Glennon, B.; Kuentz, M.; Steele, G.; Kramer, H. J. M.; Rydzak, J. W.; Chen, Z.; Morris, J.; Kjell, F.; Singh, R.; Gani, R.; Gernaey, K. V.; Louhi-Kultanen, M.; O’Reilly, J.; Sandler, N.; Antikainen, O.; Yliruusi, J.; Frohberg, P.; Ulrich, J.; Braatz, R. D.; Leyssens, T.; von Stosch, M.; Oliveira, R.; Tan, R. B. H.; Wu, H.; Khan, M.; O’Grady, D.; Pandey, A.; Westra, R.; Delle-Case, E.; Pape, D.; Angelosante, D.; Maret, Y.; Steiger, O.; Lenner, M.; AbbouOucherif, K.; Nagy, Z. K.; Litster, J. D.; Kamaraju, V. K.; Chiu, M.-S. I
DOI: 10.1021/acs.cgd.5b00720 Cryst. Growth Des. XXXX, XXX, XXX−XXX
Crystal Growth & Design
Article
(21) Andronis, V.; Zografi, G. Crystal nucleation and growth of indomethacin polymorphs from the amorphous state. J. Non-Cryst. Solids 2000, 271, 236−248. (22) Wu, J. X.; Xia, D.; van den Berg, F.; Amigo, J. M.; Rades, T.; Yang, M.; Rantanen, J. A novel image analysis methodology for online monitoring of nucleation and crystal growth during solid state phase transformations. Int. J. Pharm. 2012, 433, 60−70. (23) Larsen, P. A.; Rawlings, J. B.; Ferrier, N. J. An algorithm for analyzing noisy, in situ images of high-aspect-ratio crystals to monitor particle size distribution. Chem. Eng. Sci. 2006, 61, 5236−5248. (24) Van Eerdenbrugh, B.; Baird, J. A.; Taylor, L. S. Crystallization tendency of active pharmaceutical ingredients following rapid solvent evaporationclassification and comparison with crystallization tendency from undercooled melts. J. Pharm. Sci. 2010, 99, 3826−3838. (25) Van Eerdenbrugh, B.; Taylor, L. S. An ab initio polymer selection methodology to prevent crystallization in amorphous solid dispersions by application of crystal engineering principles. CrystEngComm 2011, 13, 6171−6178. (26) Gonzalez, R. C.; Woods, R. E. Digital image processing, 3rd ed.; Prentice Hall: Upper Saddle River, NJ, 2008. (27) Barrett, H. H., III The Radon Transform and Its Applications. In Progress in Optics; Wolf, E., Ed.; Elsevier, 1984; Vol. 21, pp 217−286. (28) Moré, J. The Levenberg-Marquardt algorithm: Implementation and theory Numerical Analysis; Watson, G., Ed.; Springer: Berlin/ Heidelberg, 1978; Vol. 630, pp 105−116. (29) Xia, D.; Ouyang, M.; Wu, J. X.; Jiang, Y.; Piao, H.; Sun, S.; Zheng, L.; Rantanen, J.; Cui, F.; Yang, M. Polymer-Mediated Antisolvent Crystallization of Nitrendipine: Monodispersed Spherical Crystals and Growth Mechanism. Pharm. Res. 2012, 29, 158−169. (30) Wu, J. X.; Yang, M.; Berg, F. v. d.; Pajander, J.; Rades, T.; Rantanen, J. Influence of solvent evaporation rate and formulation factors on solid dispersion physical stability. Eur. J. Pharm. Sci. 2011, 44, 610−620. (31) Lindfors, L.; Forssén, S.; Westergren, J.; Olsson, U. Nucleation and crystal growth in supersaturated solutions of a model drug. J. Colloid Interface Sci. 2008, 325, 404−413. (32) Taylor, L. S.; Zografi, G. Spectroscopic Characterization of Interactions Between PVP and Indomethacin in Amorphous Molecular Dispersions. Pharm. Res. 1997, 14, 1691−1698.
J
DOI: 10.1021/acs.cgd.5b00720 Cryst. Growth Des. XXXX, XXX, XXX−XXX