1.1.

Cytoskeleton

The cytoskeleton is a network of proteins that structurally and dynamically organize the cytoplasm of living cells.¹ It is composed of three major structural elements: microtubules, intermediate filaments, and microfilaments, each consisting of polymers of protein subunits.

The cytoskeleton is responsible for the maintenance of cell architecture, shape, and internal organization. For example, it enables the transport of organelles and vescicles. Microtubules and actin filaments play a role in mitosis, cell signaling, and motility.²

A well-organized cytoskeletal network exhibits, as in the top tile of Fig. 1 and in Fig. 2 , bundles of microtubules radiating from the microtubule organizing center (MTOC), located near the cell nucleus, toward the cell periphery.

Fig. 1

Outline of the experiment design. All images show microtubules of cultured rate hepatocytes obtained by immunolocalization and epifluorescence microscopy. $C$ , tile from the untreated (control) set; $T\mathit{502}h$ , tile from treatment at the highest concentration $(50\mu \mathrm{g}∕\mathrm{ml})$ for the longest time $\left(2\mathrm{h}\right)$ ; $R$ ; tile from the recovery experiment; and $T\mathit{2530}$ ; tile from an intermediate treatment. In tiles $T\mathit{502}h$ and $T\mathit{2530}$ the dark round area corresponds to the nucleus, in proximity of which the MTOC is located. Square side $\text{length}=13\mu \mathrm{m}$ .

Fig. 2

Microtubules of a control cell, tile conṯ03o. Square side $\text{length}=13\mu \mathrm{m}$ .

In general, alterations of cytoskeletal functions are related to (may be the cause or the effect of) a broad spectrum of biochemical reactions such as adenosine triphosphate (ATP) depletion, disruption of intracellular ion $\left({\mathrm{Ca}}^{2+}\right)$ homeostasis, thiol oxidation, and phosphorylation,^{3, 4, 5, 6, 7} all of which have major repercussions on cell functionality. Many substances are known to directly or indirectly interact with cytoskeletal constituents and cause damage: metals,^{8, 9} herbicides and fungicides,^{3, 10} the neurotoxin⁴ MPTP, as well as natural toxins.⁵

In morphological terms, damage to the microtubule network consists of structural disorganization, inhibition of assembly, disassembly, or depolymerization. All of these patterns can be visualized by means of suitable techniques¹¹ (e.g., Sec. 3).

In turn, morphological alterations of the cytoskeleton reflect functional disturbances of the whole cell, e.g., altered transport of very low density lipoproteins in hepatocytes,¹² neurodegeneration,^{13, 14, 15} and cell transformation.¹⁶ For this last reason the cytoskeleton has become one of the preferred targets of anticancer drugs.^{17, 18, 19, 20, 21, 22}

Therefore cytoskeletal morphology is a valuable indicator of cell functionality. In this context, morphological analysis and classification of cytoskeletal images can assist in estimating the degree of cell injury.

1.2.

Image Classification

Image classification can be either dichotomic or of “continuous” type. The former is meant to provide yes/no answers, the latter is expected to rank images according to some predefined properties. Moreover, image classification paradigms differ by the type of data processed. Raw images can be supplied as such to a black box algorithm, which is requested to sort out similarities and rank images accordingly.^{23, 24} Instead, the paradigm implemented in this paper starts with the extraction of a few quantitative morphological descriptors, also called indicators or features, believed to effectively translate structural, hence functional information. Features were extracted by the following methods: contour (Secs. 5.2, 6.1) and mass fractal analysis (Secs. 5.3, 6.2); spatial differentiation (Secs. 5.4, 6.3); and “spectrum enhancement,” a nonlinear filter in the Fourier domain (Secs. 5.5, 6.4). Descriptor selection is based on known function-morphology relations and driven by the properties of the available images (Sec. 4.1).

Classification and quantitative assessment then relies on the four standard steps:

These stages are all reflected in the experiment design (Sec. 4) and the developments described in the following, where the eventual application is the quantitative description of cytoskeletal morphology affected by damage (Sec. 9.1) or undergoing recovery (Sec. 9.2).

4.1.

Required Morphological Descriptors

Since experiments were carried out in different situations, the average fluorescence yield affecting a given image varied. Moreover, the photograph sets of a given experiment were developed and printed at different times. As a consequence, images had to be analyzed by methods exhibiting the least possible sensitivity to the average fluorescence intensity and to the absolute relation between exposure and film density.

The features to be extracted from the images for classification had to translate into quantitative terms those properties that a human expert believes to characterize cytoskeletal organization or loss thereof. In morphological terms, said features depend on image structure (leading features) and, to some extent, on texture (fine details).

The fully developed microtubules of a normal cytoskeleton are known to have indented contours, whereas an altered cytoskeleton has a smoother contour. One morphological descriptor known to quantify indentation is the contour fractal dimension, denoted by $D_{\mathrm{FP}}$ (Secs. 5.2, 6.1).

The number of microtubule bundles per unit area, i.e., the “mass density,” of a normal cytoskeleton is lower than that of an altered one. A suitable descriptor, suggested, e.g., by the morphological analysis of dendritic materials,⁴⁶ is the mass fractal dimension $D_{\mathrm{FM}}$ (Secs. 5.3, 6.2).

Microtubule bundles give rise to more or less visible edges in a gray-scale image. Global, direction-independent edge indicators are therefore necessary, such as the total variation (TV) and a suitable norm of the Laplace operator $L$ . Both descriptors are computed in the direct domain, i.e., by (discrete) spatial differentiation (Secs. 5.4, 6.3).

Finally, a set of descriptors was required that

Separation of image structure from texture is one of the fundamental tasks of image understanding, which was formalized by the Osher-Rudin paradigm (Chap. 1 of Ref. 47). In this application, structure observed in a cytoskeleton is due to organized microtubule bundles, whereas (straight or curly) isolated tubules and image noise contribute to texture. The diameter of a microtubule is $25\mathrm{nm}$ . If straight bundles are formed by, say, 5 to 20 microtubules, the image will contain structures of spatial period ranging from $250\text{to}1000\mathrm{nm}$ . The corresponding fundamental spatial frequency will range from $4\text{to}1\text{cycles}∕\mu \mathrm{m}$ . Hence the analysis of the power spectrum shall be tuned to that frequency domain.

The implementation of these requirements by means of “spectrum enhancement” is described in Secs. 5.5, 6.4. Note also that Table 1 summarizes the symbols used in this paper.

Table 1

List of symbols.

4.2.

Classification Experiment Design

The classification experiment was designed to implement all three stages, i.e., training, validation and recognition. Figure 1 provides examples of the image types and an outline and Table 2 the complete list of processed tile sets. Untreated cells, which underwent all sample preparation stages described in the previous section except exposure to Benomyl™, served as negative control: they yielded comparable images divided into the $C$ and $C_V$ tile sets for experimental purpose. The $T\mathit{502}h$ and $T\mathit{502}{h}_V$ sets were derived from treatment with Benomyl™ at experimental conditions( $50\mu \mathrm{g}∕\mathrm{ml}$ for $2\mathrm{h}$ ) classified as “extreme” according to the morphological and biochemical evidence given by Ref. 3. Cells from intermediate treatments gave rise to other tile sets, e.g., $T\mathit{2530}$ ( $25\mu \mathrm{g}∕\mathrm{ml}$ for $30\min$ ). Cells from $24\text{-}\mathrm{h}$ recovery after the $T\mathit{502}h$ treatment yielded the $R$ tile set. The $C$ and $T\mathit{502}h$ sets were used for classifier training (Sec. 7) as well as the $X$ tile (Sec. 3 and Fig. 6 in Sec. 6.4). The $C_V$ and $T\mathit{502}{h}_V$ tiles were used for classifier validation (Sec. 8). Sets $T\mathit{5015}$ , $T\mathit{5030}$ , $T\mathit{252}h$ , and $T\mathit{2530}$ were used in the recognition stage to rank the effects of intermediate concentrations and exposure times on morphology (Sec. 9.1). The analysis of set $R$ is presented in Sec. 9.2.

Table 2

Image sets forming the experiment design.

5.1.

Preprocessing and Optimal Thresholding

Given the heterogeneity of experimental conditions and of photographic processing (Sec. 4.1), independence of the absolute fluorescence intensity could be achieved by normalizing the gray-scale histogram of each tile separately. Histogram normalization is an affine transformation of $g\left(\mathbf{x}\right)$ into $f\left(\mathbf{x}\right)$ according to

The extraction of a few morphological descriptors required image thresholding, e.g., for binarization or background subtraction. Several methods exist that determine a binary threshold according to some optimality criterion. In this paper, the optimum threshold for each tile $\tau$ was determined to maximize cross-correlation ${\rho }_{g\gamma }\left(t\right)$ between the original gray scale $g$ and the black-and-white image $\gamma \left(t\right)$ binarized according to the threshold $t$ . Formally,

5.2.

Contour Fractal Analysis

The purpose of contour fractal analysis is to estimate the corresponding fractal dimension $D_{\mathrm{FP}}$ . A typical algorithm consists of the following eight stages: (1) $\tau$ thresholding (Sec. 5.1); (2) determination of the Férét diameter $F_D$ ; (3) contour tracking; (4) estimation of the perimeter $P\left(y\right)$ by different values of the yardstick $y$ ; ranging from $y=(1∕2)F_D$ to $y=1\text{pixel}$ ; (5) contour correction wherever applicable; (6) construction of the $\log \left[P\left(y\right)\right]$ versus $\log \left(y\right)$ Richardson plot; (7) determination of the morphological threshold $y_{\mathrm{MT}}$ ; and (8) linear interpolation of the plot in $1⩽y⩽y_{\mathrm{MT}}$ . Some of the computer algorithms can be found in Ref. 48.

As shown in Appendix B, Proposition 1, the value of $D_{\mathrm{FP}}$ estimated from an image binarized by the $\tau$ of Eq. 2 is invariant with respect to affine transformations of $g(\bullet )$ such as histogram normalization [Eq. 1]. As a consequence $D_{\mathrm{FP}}$ is an intrinsic property of the analyzed tile, not affected by absolute gray levels.

5.3.

Mass Fractal Analysis

Among the methods that estimate the mass fractal dimension $D_{\mathrm{FM}}$ , box counting⁴⁶ is the most straightforward. It also operates on a binarized image and relies on regression of a Richardson plot. Details are left to the implementation (Sec. 6.2).

5.4.

Direct Methods (Spatial Differentiation)

Although all methods are applied to a discrete context, a continuum setting now simplifies the description. By assuming

5.5.

“Spectrum Enhancement”

The spectrum enhancement algorithm, suggested by the experimental results from a coherent optical image processor many years ago⁴⁹ and specifically developed^{41, 50, 51, 52, 53} for image classification in the past $3\mathrm{yr}$ , meets all the specifications listed in Sec. 4.1.

Let $Q\Omega$ denote a square of side-length $L$ and consider an image, i.e., a function $Qg\left(\mathbf{x}\right)$ , $\mathbf{x}\equiv \{x_1,x_2\}∊Q\Omega$ , which is continuous on the surface $\mathcal{T}$ of the torus obtained by glueing the opposite sides of $Q\Omega$ together. One way of obtaining such a $Qg(\bullet )$ from an image $g(\bullet )$ defined in a square tile $\Omega$ of side-length $L∕2$ is the application of the twofold $Q(\bullet )=\text{flop}\left[\text{flip}(\bullet )\right]$ reflection. Next let $Q\Omega$ be discretized by a square grid of step-length $\mathcal{l}$ . Let $\mathbf{u}\equiv \{u_1,u_2\}$ be the spatial frequency vector. Then the discrete Fourier transform $G\left(\mathbf{u}\right)$ of $Qg(\bullet )$ is defined in the square

Represent $\mathbf{u}$ in polar coordinates $\mathbf{u}\equiv \{u,\theta \}$ , where $u=\mid \mathbf{u}\mid$ is the wave number, and $\theta$ is the polar angle such that $0⩽\theta ⩽2\pi$ . Denote by ${\mid G\left(\mathbf{u}\right)\mid }^2$ the power spectral density. Due to the symmetry of $Qg(\bullet )$ , the main axes of inertia of ${\mid G\left(\mathbf{u}\right)\mid }^2$ (those that diagonalize the inertia tensor) coincide with ${\mathbf{u}}_1$ and ${\mathbf{u}}_2$ . Let $\Theta$ denote an arc symmetric with respect to either axis. Then the normalized, arc-averaged spectral density profile is the function $s(\bullet )$ of $u$ alone defined in $0⩽u⩽u_{\max }$ (cycles/image) according to

Let $m\left(u\right)$ be a model spectral density. For example, choose

Among the properties of the enhanced spectrum $h\left[u\right]$ the following two are worth mentioning.

Proposition 2. In the interval $0⩽u⩽u_{\max }$ , $h\left(u\right)$ is invariant with respect to scaling transformations such that $\beta =0$ in Eq. 1, applied to the gray levels of the input image.

This is a straightforward consequence of normalization [Eq. 10]. As a result, the enhanced spectrum reflects the intrinsic properties of the analyzed tile.

The second property is a relation between $h\left(\mathbf{u}\right)$ and spatial differentiation of $Qg\left(\mathbf{x}\right)$ .

Proposition 3 (in words). When $p∕2$ is an integer, the enhanced spectrum corresponds to evaluating all spatial derivatives of order $p∕2$ of the image, taking their Fourier transforms, forming a linear combination of the squared moduli, adding $\delta \left(u\right)$ , taking the logarithm and averaging over $\Theta$ . The formal statement and other remarks are provided in Appendix C.

From a practical point of view, conversion of ${\mid G\left(\mathbf{u}\right)\mid }^2$ to the log scale [Eq. 10] amplifies the contribution of structured but weakly emitting microtubules. Averaging over $\Theta$ makes $s\left(u\right)$ independent of bundle direction. Comparison between the asymptotics of the given power spectrum and a model of the $m^{\left(p\right)}(\bullet )$ class has been implied by previous studies involving visual processing.⁵⁴ The emphasis herewith is on the low-frequency behavior of $h(\bullet )$ for example in $1⩽u⩽(1∕4)u_{\max }$ , aimed at enhancing image structure and assessing the relative weight of texture. Given the scale factor of Eq. 9, the approximate spatial frequency range to which microtubule bundles shall contribute most in terms of power spectral density is

6.1.

Contour Fractal Dimension $D_{\mathrm{FP}}$

Contour correction, step 5 of Sec. 5.2, consisted of subtracting from $P\left(y\right)$ the $y$ independent length of the straight segments along the boundary of the square before building the Richardson plot. Otherwise $D_{\mathrm{FP}}$ would have been underestimated. The threshold $y_{\mathrm{MT}}$ (step 7) was determined once for all from the Richardson plot of a contour of known dimension. This is a yardstick value that separates the structural part (large $y$ , highly oscillatory graph) of the plot from the textural one (small $y$ , smoother graph). The Richardson plots yielded by the tiles of Figs. 2 and 3 are shown by Fig. 4 . Contour correction was applied to the former plot. The latter tile, instead, which showed a whole cytoskeleton, did not need any correction. The estimated values of $D_{\mathrm{FP}}$ for $C$ tiles, including that of Fig. 2, were affected by the deep contour indentations due to filaments. The smoother contours of $T\mathit{502}h$ tiles (Fig. 3) yielded lower values, as expected.

Fig. 3

Microtubules of a treated cell, tile be502̱ 04. Square side $\text{length}=13\mu \mathrm{m}$ .

Fig. 4

Richardson plots related to contour fractal analysis. All length values in units of the Férét diameter. For plot conṯ03o, derived from Fig. 2, estimated $D_{\mathrm{FP}}=1.54$ after discounting for straight edges, and for plot be502̱04 derived from Fig. 3, estimated $D_{\mathrm{FP}}=1.11$ . No straight-edge correction was required.

6.2.

Mass Fractal Dimension $D_{\mathrm{FM}}$

The $D_{\mathrm{FM}}$ was estimated by box counting⁴⁶ on each previously normalized and $\tau$ binarized tile. The software used was Benoit™ 1.3 of TruSoft Inc.⁵⁶ Typically, the box scaling ratio was 1.3 and the grid orientation angle was stepped by $15\mathrm{deg}$ .

The values of tile set-averaged $D_{\mathrm{FP}}$ and $D_{\mathrm{FM}}$ and their standard deviations are shown in Table 3 .

Table 3

Set averages of descriptors used by the implemented classifier.

6.3.

Total Variation ${\mathrm{TV}}_{\eta }$ and Laplacian Norm $L_{\eta }$

The relation between $\tau$ and $\eta$ mentioned at the end of Sec. 5.4 was

6.4.

Descriptors of the Enhanced Spectrum

From here on the wave number is assumed to be discrete and span the interval $0⩽u⩽u_{\max }$ . The enhanced spectrum can be regarded as the signature of a given tile, and as such, the graph $\{u,h\left(u\right)\}$ is an array of first-level morphological descriptors. The decision was made to extract fewer second-level descriptors according to some criterion. For any exponent $p,h\left(u\right)$ exhibits a trend over which jumps and oscillations are superimposed. In fact, the latter are amplified with respect to those of $s(\bullet )$ . For this reason the value $h\left(0\right)=0$ was reassigned according to $h\left(0\right)\leftarrow h\left(1\right)$ and $h(\bullet )$ was interpolated in $0⩽u⩽u_{\max }$ by a polynomial $q(\bullet ;d)$ of suitable degree $d$ . Four parameters eventually affected $q(\bullet ;d)$ : the axis of inertia ( $I=1$ : major axis, $I=2$ : minor axis) about which integration was carried out [Eq. 10], the arc $\Theta$ , the model exponent $p$ , and the degree $d$ .

Plots obtained from some tiles used in classifier training with $I=1$ , $\mid \theta \mid ⩽46\mathrm{deg}$ , $p=1.8$ , and $d=9$ are shown in Fig. 5 . In agreement with the estimate of Eq. 13, enhanced spectra from the $C$ set generally exhibited a local maximum in $0<u<63\text{cycles}$ /image as shown, e.g., by conṯ3b (Fig. 5). One exception was tile conṯ1a. None of the $T\mathit{502}h$ tiles (plots be502̱3a and be502̱1a of Fig. 5) had such maximum; in fact, $q^{\prime }(\bullet )<0$ in the same frequency range. The plot conṯx of the perinuclear tile $X$ seemed to deviate from all other plots of either set. Tile $X$ itself is shown in Fig. 6 .

Fig. 5

Interpolated enhanced spectra of representative tiles used in classifier training. Integration is about the major axis of inertia over an arc of $\pm 46\mathrm{deg}$ ; model exponent $p=1.8$ ; degree of polynomial $d=9$ . Plots conṯ1a and conṯ3b represent the $C$ set, plots be502̱1a and be502̱3a represent the $T\mathit{502}h$ set; and plot conṯx is the perinuclear tile (see Fig. 6). Local maxima in the spatial frequency range $20<u<100\text{cycles}$ /image correspond to structured microtubule bundles of a normal cytoskeleton.

The selection of second-level descriptors was a part of classifier training (Sec. 7) and initially involved the following five quantities:

The values of $u_H$ and $u_C$ were also determined in the training stage. As one can infer from the last three columns of Table 3, the first three descriptors, $A$ , $q^{\prime }\left(0\right)$ , and $c$ were eventually selected for image classification.

8.1.

Sensitivity to Descriptors

In an early training attempt,⁴⁴ the set

Fig. 8

Classification of the $C$ set (filled boxes), the $T\mathit{502}h$ set (filled triangles), and $X\left({}^{*}\right)$ by descriptor set ${\mathcal{D}}_8$ . Modulo reflections of both axes, the pattern is similar to that of Fig. 7a, in spite of the different number ( $N=8$ instead of 7) and list of descriptors [Eq. 19 instead of Eq. 16]. Note the different scales of ${\mathbf{z}}_1$ and ${\mathbf{z}}_2$ . The first principal component explains 41% of the sample variance, the first two together explain 75%, and the first three 87%.

The sensitivity of the classification to changes in the descriptor sets (from ${\mathcal{D}}_8$ to ${\mathcal{D}}_7$ ) can be now quantified as follows. The set ${\mathcal{D}}_7$ differs from ${\mathcal{D}}_8$ by three descriptors out of eight. The situation in $\{{\mathbf{z}}_1;{\mathbf{z}}_2\}$ can be summarized by the following quotients: $[\mathrm{dist}(C,T\mathit{502}h)∕d_{\mathrm{CT}}]$ , $[\mathrm{diam}\left(C\right)∕d_{\mathrm{CT}}]$ , and $[\mathrm{diam}\left(T\mathit{502}h\right)∕d_{\mathrm{CT}}]$ , where dist(∙,∙) is the distance between two clusters (minimum distance between two points belonging to either cluster), diam(∙) is the diameter of a cluster, and $d_{\mathrm{CT}}$ serves as a yardstick in the $\{{\mathbf{z}}_1;{\mathbf{z}}_2\}$ plane. The values are given by Table 4 . Ideally $[\mathrm{dist}(C,T\mathit{502}h)∕d_{\mathrm{CT}}]$ would be close to 1, and the other two will be much smaller than 1. Neither ${\mathcal{D}}_8$ nor ${\mathcal{D}}_7$ provide the “ideal” classification. However, the focus is on classifier sensitivity: therefore, a comparison between ${\mathcal{D}}_8$ and ${\mathcal{D}}_7$ has to be made in relative terms. The normalized distance between clusters $[\mathrm{dist}(C,T\mathit{502}h)∕d_{\mathrm{CT}}]$ is almost insensitive to the descriptor set: this means sets $C$ and $T\mathit{502}h$ are as “separable” by ${\mathcal{D}}_7$ , as they were by ${\mathcal{D}}_8$ . Instead, the normalized diameters $[\mathrm{diam}(\bullet )∕d_{\mathrm{CT}}]$ of both clusters $C$ and $T\mathit{502}h$ are more sensitive to the replacement of ${\mathcal{D}}_8$ by ${\mathcal{D}}_7$ , because they increase by 26 and 13%, respectively; ${\mathcal{D}}_7$ makes both training sets look “more heterogeneous” than they were under ${\mathcal{D}}_8$ . From this numerical experiment, one can conclude that the classifier is relatively “stable.” Further evidence of classifier “stability” is provided by the values of $V_k$ , $k=1,2,3$ , given in the caption of Fig. 8; the first three principal components derived from either ${\mathcal{D}}_8$ or ${\mathcal{D}}_7$ carry the same amount of information.

Table 4

Indicators of classifier sensitivity with respect to descriptor sets.

8.2.

Internal Validation

Internal validation consisted of submitting to the trained classifier some new tiles belonging to either class, then projecting the corresponding feature vectors onto $\{{\mathbf{z}}_1;{\mathbf{z}}_2\}$ by means of $\mathsf{T}$ of Eq. 17 and finally computing sensitivities and specificities.

In detail, the 8 $C_V$ and 8 $T\mathit{502}{h}_V$ tiles defined in Sec. 4 were selected from boundary areas of the available images. Unlike those used in training, some $C_V$ and $T\mathit{502}{h}_V$ tiles could not be easily classified visually. For example Figs. 9a and 9b exhibit “borderline” morphologies; at the top right of Fig. 9a, near the cell nucleus, is a dense mass of tubulin, which looks disorganized, whereas the boundary tubules of Fig. 9b are sufficiently structured as in a control cytoskeleton.

Fig. 9

(a) Tile of the $C_V$ validation set, where the dense, poorly structured microtubules on the upper right, close to the perinuclear region, look treated. Misclassification is expected, as shown in Fig. 10b. Square side $\text{length}=13\mu \mathrm{m}$ . (b) Tile of the $T\mathit{502}{h}_V$ validation set, where microtubules on the boundary, although randomly oriented, exhibit sufficient structure as in a control cytoskeleton. Misclassification is expected, as shown in Fig. 10b. Square side $\text{length}=13\mu \mathrm{m}$ .

Fig. 10

(a) Centroids of the training sets $\{C,T\mathit{502}h\}$ (gray) and of the validation sets $\{C_V,T\mathit{502}{h}_V\}$ (black). The displacement of the $T\mathit{502}{h}_V$ centroid with respect to the $T\mathit{502}h$ centroid is relatively small. (b) Validation result. Tiles of the $C_V$ set are represented by filled boxes, those of the $T\mathit{502}{h}_V$ set by filled triangles. Training tiles (Sec. 7) are shown by gray boxes and triangles respectively. Not all $C_V$ tiles are on the left half-plane as well as not all $T\mathit{502}{h}_V$ tiles are on the right one.

Descriptors were extracted by the procedure outlined in Sec. 7. By assumption, the training tiles represented their respective sets, therefore normalization consisted of subtracting the ${\mu }_n$ and dividing by the ${\sigma }_n$ computed in the preceding. The normalized feature vectors ${\mathbf{y}}_m$ , $1⩽m⩽M=16$ , were determined for each tile of the validation sets. The ${\mathbf{y}}_m\mathrm{s}$ were projected onto $\{{\mathbf{z}}_1;{\mathbf{z}}_2\}$ according to

In terms of individual tiles the result is shown by Fig. 10b. References to a few images are also displayed. The diameters of both $C_V$ and $T\mathit{502}{h}_V$ clusters are larger than those of $C$ and $T\mathit{502}h$ . Two $C_V$ tiles, one of which is Fig. 9a, are positioned closer to the $T\mathit{502}h$ set and one $T\mathit{502}{h}_V$ tile [Fig. 9b] is closer to the $C$ set. The $z_1$ coordinates of these three tiles [Fig. 10b], which deviate from those of the remaining tiles, can be explained by means of the $z_1$ loadings of the morphological indicators [Fig. 7b]. For example, the tile of Fig. 9b has $D_{\mathrm{FP}}=1.475$ , the highest value in the $T\mathit{502}{h}_V$ set, comparable to the $C$ and $C_V$ set averages (Table 3). Since $D_{\mathrm{FP}}$ is negatively correlated to $z_1$ , this property suffices to position the tile in the left half-plane. About the $C_V$ tile of Fig. 9a, its enhanced spectrum, like that of a $T\mathit{502}h$ tile, has no local maximum. The presence or absence of the latter affects all three indicators $A$ , $q^{\prime }\left(0\right)$ , and $c$ , which negatively correlate to $z_1$ . Therefore, this tile shall lie away from the centroid of the $C$ -(training) set and possibly in the right half-plane.

If $\{{\mathbf{z}}_1;{\mathbf{z}}_2\}$ is regarded as a feature plane, then inspection of Fig. 10b suggests the training sample is separable; there is a straight line, e.g., $z_1=-0.036$ , that divides the $C$ and $T\mathit{502}h$ sets. Each tile of the validation set could be classified accordingly. The confusion matrix was formed and the sensitivity and specificity values of Table 5 were obtained.

Table 5

Sensitivity and specificity of the classifier in validation mode.

8.3.

Discriminant Function Analysis

Classification is said to be supervised whenever the class to which a feature vector belongs is known beforehand. Given $K$ classes, discriminant analysis (DA for short; e.g., Chap. 3 of Ref. 57) implements supervised classification by looking for $K-1$ affine combinations, called discriminant functions, of the $N$ features (the morphological descriptors), which maximize separation between classes and minimize within-class scatter. Three classes of belonging were defined: (1) the singleton $X$ , (2) $\widehat{C}$ made of 8 $C$ and 8 $C_V$ tiles, and (3) $\widehat{T}$ made of 8 $T\mathit{502}h$ and 8 $T\mathit{502}{h}_V$ tiles, hence $M=33$ . The corresponding $N\times M$ raw feature matrix was normalized rowwise as in Sec. 7 and DA carried out both by SPSS™ and by the LDA module of Q-Parvus.⁶⁰ The $3\times 3$ classification (or “confusion”) matrix $\mathsf{K}$ was constructed. Since 14 $\widehat{C}$ tiles were assigned to the $\widehat{C}$ class and 2 to the $\widehat{T}$ class, whereas $X$ and all $\widehat{T}$ tiles were correctly classified, the relevant sensitivity (normalized row-wise sums of entries of $\mathsf{K}$ ) and specificity (normalized column-wise sums of entries of $\mathsf{K}$ ) of Table 6 resulted. A cutout of the $\{{\mathbf{D}}_1;{\mathbf{D}}_2\}$ discriminant functions plane is shown by Fig. 11 . For the purpose of classifier training and validation, DA answers the question of how suitable are the given morphological descriptors to discriminate between the given $K$ classes.

Fig. 11

Discriminant analysis where the $\widehat{C}$ and $\widehat{T}$ sets are made of training and validation tiles together; $\widehat{C}$ (boxes) and $\widehat{T}$ (triangles) tiles are represented in a plane where the discriminant functions $\{{\mathbf{D}}_1;{\mathbf{D}}_2\}$ form the reference frame. The $X$ tile is located further away from the origin.

Table 6

Sensitivity and specificity obtained from discriminant analysis.

Sensitivities and specificities shall all be greater than 50%, regardless of $K$ . In the ideal case, they are all equal to 100%, which corresponds to a diagonal $\mathsf{K}$ . From the entries of Table 6 one deduced that information stored in the chosen descriptor set, ${\mathcal{D}}_7$ , could adequately and reliably tell $\widehat{C}$ from $\widehat{T}$ tiles, a necessary condition for going over to the next, more challenging tasks: recognizing and ranking tiles from other experiments.

9.1.

Effects of Different Treatments

Projection onto $\{{\mathbf{z}}_1;{\mathbf{z}}_2\}$ scattered the $T\mathit{5015}$ tiles almost evenly in the rectangle $\{-2.7⩽z_1⩽2.8,-0.9⩽z_2⩽3.0\}$ , as shown by Fig. 12 . Their centroid lay at ${\Gamma }_{T\mathit{5015}}\equiv \{-0.082,0.924\}$ . Tiles of the $T\mathit{5030}$ set (Fig. 9 of Ref. 45) were mostly aggregated in the $T\mathit{502}h$ area. Their centroid lay at ${\Gamma }_{T\mathit{5030}}\equiv \{1.529,-0.362\}$ .

Fig. 12

$T\mathit{5015}$ tiles (filled triangles) in the $\{{\mathbf{z}}_1;{\mathbf{z}}_2\}$ plane. Tiles derived from exposure to $50\mu \mathrm{g}∕\mathrm{ml}$ for $15\min$ . Training $C$ (gray boxes) and $T\mathit{502}h$ tiles (gray triangles) are shown for reference.

Classification of tiles derived from exposure at $25\mu \mathrm{g}∕\mathrm{ml}$ is represented by Fig. 13 , which pertains to the $T\mathit{2530}$ set. Comments are included in the caption. Results for the $T\mathit{252}h$ set were reported in Fig. 10 of Ref. 45. The centroid coordinates were ${\Gamma }_{T\mathit{252}h}\equiv \{1.714,-0.215\}$ and ${\Gamma }_{T\mathit{2530}}\equiv \{1.028,-0.721\}$ , respectively.

Fig. 13

Tiles $T\mathit{2530}$ (filled disks) in the $\{{\mathbf{z}}_1;{\mathbf{z}}_2\}$ plane. Tiles derived from exposure to $25\mu \mathrm{g}∕\mathrm{ml}$ for $30\min$ . Training $C$ (gray boxes) and $T\mathit{502}h$ tiles (gray triangles) are shown for reference. As compared to $T\mathit{252}h$ tiles, $T\mathit{2530}$ tiles are more scattered. Some are closer to the $C$ half-plane. The tile placed at the extreme left exhibited the typical microtubule structure of a control cytoskeleton.

9.2.

Quantitative Estimate of Structural Recovery after Treatment

Undisputable experimental evidence³ showed that cytoskeletal structures are capable of recovering after exposure to toxic substances, provided that dose does not exceed some threshold. Visual inspection of $R$ images (Sec. 3) suggested that cytoskeletal damage had to be reversible, at least in part. The first quantitative result about structural recovery was obtained by the ${\mathcal{D}}_8$ classifier (Ref. 44 and Sec. 8.2). Herewith, the ${\mathcal{D}}_7$ classifier (Sec. 7) was applied to the $R$ tiles and yielded the result of Fig. 14 . Only 3 out of 12 tiles are located in the $T\mathit{502}h$ half-plane and the remainder in the $C$ half-plane.

Fig. 14

Tiles $R$ (filled diamonds) in the $\{{\mathbf{z}}_1;{\mathbf{z}}_2\}$ plane. This result provides evidence of 0.73 $(+0.18,-0.07)$ structural recovery $24\mathrm{h}$ after extreme treatment ( $50\mu \mathrm{g}∕\mathrm{ml}$ for $2\mathrm{h}$ ). Tiles $C$ and $T\mathit{502}h$ are shown for reference.

From the $z_1$ coordinates of centroids, $z_{1,T\mathit{502}h}$ , $z_{1,R}$ , and $z_{1,C}$ one can form the quotient

9.3.

Ranking by Centroids

To summarize the results, the centroids of all sets, except those used in validation, are displayed by Fig. 15 . According to the $z_1$ coordinate alone, the sets form the sequence

Fig. 15

Set centroids, where ${\mathfrak{R}}_{\text{centroid}}=\{-2⩽z_1⩽2,-1⩽z_2⩽1.5\}=\text{rectangle}$ containing the centroids and ${\mathfrak{R}}_{\text{recognition}}=\{-5⩽z_1⩽3.5,-2⩽z_2⩽3\}=\text{rectangle}$ circumscribed to all classified tiles. Comments are provided in Sec. 10.3. The $X$ tile is not shown.

10.1.

Image Formation, Capture, and Analysis

All images were obtained from an ordinary epifluorescence microscope, acquired from photographic film developed and printed after each experiment. Therefore,

On one hand, the wide dynamic range of film obviously was an advantage because it made weakly emitting microtubules visible. On the other hand, parameter variability would have prevented classification based on absolute fluorescence values. The threshold of Eq. 2 (which implies Proposition 1, Appendix B) and the spectrum enhancement algorithm of Secs. 5.5, 6.4 (to which Proposition 2 applies) were selected as countermeasures for 1 and 2.

The opportunity of replacing ordinary by confocal microscopy in the morphological analysis of the cytoskeleton is still a matter of discussion; in Ref. 38, confocal images were successfully processed; however, there are situations⁶¹ where ordinary microscopy has been more informative.

In principle the whole procedure, from analysis to classification, can be applied to images acquired by different microscopes and sensors.

10.2.

Meaning of Some Morphological Descriptors

The estimated values of the fractal dimensions $D_{\mathrm{FP}}$ and $D_{\mathrm{FM}}$ agree with expectations (Sec. 4.1) and are easily interpreted accordingly. Fully developed microtubules of the control set do have indented contours (higher $D_{\mathrm{FP}}$ ) and lower mass density $\left(D_{\mathrm{FM}}\right)$ . The loss of organization caused by treatment (sets $T\mathit{502}h$ , $T\mathit{252}h$ , and $R$ ) yields a more compact cytoskeleton (higher $D_{\mathrm{FM}}$ ) with a smoother contour (lower $D_{\mathrm{FP}}$ ).

The two descriptors from direct methods (Secs. 5.4, 6.3) are linked to visual properties of images. For example, the well organized microtubules of $C$ cells (Fig. 2) suggest a higher value of total variation ${\mathrm{TV}}_{\eta }$ than that of $T\mathit{502}h$ cells (Fig. 3), which is in agreement with Table 3. The norm of the Laplacian, $L_{\eta }$ , is uncorrelated to $z_1$ [Fig. 7b] and uniquely characterizes $X$ : microtubules around the MTOC, hence in tile $X$ , appear entangled under an ordinary, nonconfocal microscope, and as such do not exhibit any symmetry. This justifies the very large $L_{\eta }$ of $X$ and the anticorrelation of $L_{\eta }$ to $z_2$ . If classes are ranked by increasing ${\mathrm{TV}}_{\eta }$ one obtains the following sequence (Table 3)

In spectrum enhancement (Secs. 5.5, 6.4) the subtraction of $m^{\left(p\right)}$ [Eq. 11] was the key feature of the algorithm. Namely, the subtraction of $m^{\left(p\right)}(\bullet )$ from $s(\bullet )$ [Eq. 12] represents the deviation of $s(\bullet )$ from the proposed asymptotics. Moreover, polynomial interpolation reduced the sensitivity of morphological indicators to the oscillations of $h(\bullet )$ . Separation of structure from texture is context dependent. In this application the values of either $h(\bullet )$ or $q(\bullet )$ in the interval $0⩽u⩽63\text{cycles}$ /image (approximately $2500\text{cycles}∕\mathrm{mm}$ in the case of Figs. 2 and 3) represented image structure, whereas those in $63⩽u⩽u_{\max }∕2$ corresponded to those textural details, which contributed to classification. In particular, local maxima of $q(\bullet )$ in $0⩽u⩽63\text{cycles}$ /image were originated by high contrast bundles of microtubules, typical of control cells, not observed in damaged cytoskeletal proteins (Fig. 3). Since the selected degree of $q(\bullet )$ is relatively low $(d=9)$ , it is obvious that a local maximum to the right of the origin correlates with both $A>0$ and $q^{\prime }\left(0\right)>0$ . If said maximum is peaked enough, then the $c$ of Sec. 6.4 is positive (lines $C$ and $C_V$ of Table 3). In the sets $T\mathit{502}h$ and $T\mathit{502}{h}_V$ all three descriptors have opposite signs. The simultaneous occurrence of these properties explains why $A$ , and $q^{\prime }\left(0\right)$ have very close factor loadings [Fig. 7b].

The angle-averaged spectra shown in Refs. 36 (logarithmic intensity scale) and 37 (linear scale) can be compared to $s(\bullet )$ of Eq. 10. They do exhibit local maxima even without enhancement.

The remarks in this section should help in justifying the classification result of Fig. 7a.

10.3.

Overall Classifier Performance

The statistical significance of the described classification results is supported by the following arguments and figures. Section 8 provided the necessary evidence of classifier “stability” and validity for continuing on to applications. Recognition of a relatively large set of new tiles (67) in Sec. 9 were essentially an extended validation test of a classifier trained by a total of $M=17$ tiles: namely, the $R$ , $T\mathit{5015}$ , $T\mathit{2530}$ , $T\mathit{5030}$ , and $T\mathit{252}h$ centroids (Fig. 15) all belong to the rectangle ${\mathfrak{R}}_{\text{centroid}}=\{-2⩽z_1⩽2,-1⩽z_2⩽1.5\}$ , i.e., they fall within the strip $\{-2⩽z_1⩽2\}$ defined by the $C$ and $T\mathit{502}h$ centroids after training. Moreover, none of the recognized tiles falls outside the rectangle ${\mathfrak{R}}_{\text{recognition}}=\{-5⩽z_1⩽3.5,-2⩽z_2⩽3\}$ and 46 out of 67, i.e., 68.6% fall within ${\mathfrak{R}}_{\text{training}}$ of Eq. 18. These results should not be taken for granted and prove that training based on a set of 17 tiles has captured enough of cytoskeletal morphology.

The correspondence between the gross a priori threefold distinction

10.4.

Innovative Features

All image analysis methods quoted in Sec. 2 represented significant advances toward automated classification. However, each of them relied on only one type of morphological descriptor. Instead, this paper implemented the fusion of heterogeneous morphological descriptors: fractal analysis, direct methods, and spectrum enhancement together have led to the fingerprint of an image.

Moreover, spectrum enhancement is a relatively new algorithm. It was introduced by the first author and so far it has been successfully applied to classify scattering patterns⁵⁰ and electron microscope images of other materials such as ceramic nanoaggregates^{51, 52} and tire tread particles.⁵³ Its first application to images of cytoskeleta was presented in Ref. 41.

The third innovative feature of this paper is the application of the trained classifier to rank structural damage and quantify recovery.

10.5.

Possible Developments

Although multivariate statistics is a standard tool in image classification and pattern recognition, feature extraction remains task specific and deserves additional effort. Therefore, some improvements could be brought in to increase classifier throughput and to speed up the learning phase:

On the applications side, the analysis of image sets coming from different experiments on different cell lines would outline the “operating region” of the whole approach.