1 Introduction

The extracellular matrix (ECM) provides a biophysical and biochemical environment for cell mechanical behavior. In turn, cellular interactions with the ECM resulting from adhesive, proteolytic and migratory activity govern continuous matrix reorganization. One such example occurs in eutrophic inward remodeling of small arteries [34]. Here, the existing collagen matrix is rearranged around a smaller diameter. Such inward remodeling is a hallmark of many hypertensive disorders and has been shown to have predictive value for cardiovascular events [9].

In vitro setups of cell-seeded matrix scaffolds allow studying the specific components that are active in vascular remodeling [22]. Especially, matrix reorganization can be monitored in the presence of specific cells, matrix elements, and blockers or markers of remodeling enzymes such as the matrix metalloproteases and transglutaminases [1, 2]. Several approaches using multi-cellular preparations have been developed. Thus, free-floating collagen gels can be seeded with contractile cells, and the remodeling is monitored from the dimensions of these gels [8, 13, 16, 20]. In the so-called culture force monitor, collagen gels mixed with cells are allowed to polymerize between a fixed plate and a force transducer. Contractile force development is measured during subsequent culture of this preparation [12, 14, 22, 36]. However, interpretation of data derived from these multi-cellular approaches is troubled by variability in cellular mechanical activity, lack of information on cell density, and synergistic contractile effects [30]. Therefore, microscopic observations on matrix reorganization by single cells may provide more detailed and fundamental information on the mechanisms of remodeling.

Several groups have studied the effect of single cell activity on local deformation of matrices [5, 28]. The most frequently used substrates are silicone rubber [3, 6, 18, 21] and polyacrylamide sheets. Local deformation is then monitored by tracing the texture, sometimes facilitated by using fluorescent beads or imprinted micropatterns. Although these inert scaffolds provide a geometrical framework to study processes like cell locomotion, they do not resemble the physiological ECM. This has been solved in some cases by coating the polyacrylamide membrane with a thin layer of collagen [11, 17, 26, 27, 37]. However, remaining issues are the orientation of the collagen fibrils in the coating as compared to native collagen, and the inability of the cells to chemically modify (e.g. by cross-linking) the artificial sheets.

A possible way to overcome these issues is the use of collagen-based matrices. Cell-induced deformation patterns of these matrices and the degree of reversibility of such deformation provide information on physical remodeling and the underlying biochemical processes. Deformations have been assessed by manually tracing beads or landmarks in consecutive images [6, 28, 30]. However, so far, only a few algorithms were developed for automatic detection of substrate deformations. In particular, very few of such algorithms allow tracing in the absence of embedded beads [23, 33, 35], and more stable algorithms are required.

While single cell observations provide a more fundamental insight into matrix remodeling than macroscopic studies, a concern is efficiency. Time-lapsed video microscopy of single cells forms a labor-intensive and time-consuming experiment. Meaningful experiments require the comparison of large data sets, including different cell types, various matrix compositions, or the use of biochemical or molecular interventions. A fully automated, motorized microscopy setup is required that scans series of individual cells and their surrounding matrix. In addition, each individual cell should automatically be kept in focus during remodeling experiments, which take up 24 h or more. Finally, stable algorithms are needed for the automated analysis of geometrical reorganization with minimal user input.

The aim of this study is to develop an automated technique that can be employed to obtain and analyze local matrix remodeling by individual cells. The system that we present allows for monitoring of ∼75 cells in parallel, using time-lapse video microscopy and computer-controlled stage positioning. In addition, we present and evaluate a new algorithm for automated detection of collagen matrix deformation around these cells.

2 Methods

2.1 Cell culturing and collagen matrix preparation

Smooth muscle cells, obtained from mesenteric small arteries, were cultured in Leibovitz medium with 10% (v/v) heat-inactivated fetal calf serum. Cells from passages three to nine were used in experiments.

Matrix constructs were produced from calf skin collagen (MP Biomedicals) at a concentration of 1 mg/ml; pH was buffered by HEPES, and a mix of antibiotics (PSF and ciproxin) was added. Immediately after preparation at 4°C, the collagen mixture was poured into a 3.8 cm2 culture well and a 1.5-h polymerization period at 37°C was allowed. Then, SMCs were seeded at a concentration of about 1 cell per mm2 in the presence of 1 ml serum-free Leibovitz medium. The cells were maintained in an incubation chamber that was set to a temperature of 37°C throughout the experimental procedure. After a stabilization period of about 1 h, cell–matrix interactions were monitored by microscopic imaging for a period of 24 h of spontaneous cell contraction.

2.2 Automated microscopic imaging

In order to enable simultaneous monitoring of cell-induced matrix remodeling at multiple locations, microscopy was combined with a motorized stage. Individual cells and their surrounding matrix were studied by phase-contrast microscopy (Olympus IMT-2 with 10× objective and 2.5× projection lens). Images were captured by a Qimaging Retiga SRV camera. The calibration factor for these images (1,392 × 1,040 pixels) was 0.88 μm/pixel. The microscopic field of view was set by a motorized stage, controlled by custom written software (Matlab 7.0 with Image Acquisition Toolbox 2.0). After manually determining and storing a set of x, y, z-coordinates for about 60–90 appropriate cells, these positions were tracked through time at a 15-min interval and images were captured by an automated procedure.

During the time-lapsed image acquisition, samples were kept in focus by means of implementation into the acquisition software of one of the general auto-focus algorithms. Image contrast is optimal when a histogram of intensity values shows a broad distribution over all bins. This characteristic feature can be approximated by the standard deviation of pixel intensity values. For each image, contrast was enhanced by histogram equalization, and standard deviation was calculated. A normalized focus index (FI) was defined by dividing the standard deviation of the original image (SDoriginal) by the standard deviation of the contrast enhanced image (SDenhanced):

$$ {\text{FI}} = \frac{{{\text{SD}}_{{{\text{original}}}} }} {{{\text{SD}}_{{{\text{enhanced}}}} }} $$

For each cell, a series of images was captured at five different heights. The image was then defined to be in focus at the vertical level of maximal normalized focus. This height (z) was used as the central level when capturing the z-series in the next time step, thus allowing gradual vertical shift during the time-lapsed image acquisition. Finally, a stack of time-lapsed images was constructed for each x, y-position and analyzed off-line.

2.3 Gel dynamics analysis

Cell–matrix interactions were quantified offline using in-house designed, automated image analysis software (Matlab 7.0 with Image Processing Toolbox 4.2). Matrix reorganization was assessed by calculation of the displacement field around a cell. This was achieved by performing a cross-correlation between each two successive images in an image stack; resolution of the displacement field was refined by correlation of subimages of decreasing size. This procedure is explained below and illustrated in Fig. 1.

Fig. 1
figure 1

Graphical representation of parameters used in matrix deformation analysis. The image shows a single SMC in the center, surrounded with a collagen matrix of relatively smooth texture. The white box (768 × 768 pixels) indicates the area for the most coarse correlation analysis. This correlation between two successive images was applied over an area as indicated by the yellow box (960 × 960 pixels). The green boxes show the more refined cross-correlation windows (stages 2–5). Table 1 indicates the area expansion used for these correlation analyses

Full size image

First, gross displacement was defined at the point of maximal correlation between two parent images I 1 (at t = t 0) and I 2 (at t = t 0 + Δt), using a correlation threshold of 0.5 and a maximum tested displacement of 11 pixels. Then, I 1 was decomposed into four equal-sized square subimages, for which another cross-correlation was performed. The position of the correlation subwindow of I 2 was refined using the displacement as calculated for the parent image, with a safety measure to prevent I 2 crossing the borders of I 1. This procedure of accurately positioning a correlation window not only reduced the possibility of an accidental cross-correlation match between a subimage of I 1 and any random subimage of I 2, but also drastically reduced calculation time. Resolution of the displacement field was refined to the fifth decomposition stage in our experiments (see Fig. 1), using maximum tested displacements (n cc) as indicated in Table 1. When insufficient correlation (r < 0.5) was found in any small area, displacement field from the next larger images were used. The calculated displacements were then assigned to the corresponding subimage centers. Subsequently, a continuous displacement field was obtained by bicubic interpolation. Finally, matrix compaction for each stack was quantified by tracking circular areas centered around the cell (see Fig. 5).

Table 1 Settings for matrix compaction analysis by decomposition cross-correlation
Full size table

2.4 Validation

The method described above (decomposition CC) was validated on several test series against a straightforward cross-correlation analysis (direct CC), with settings according to decomposition stage 1.

The first test case consisted of an image of a collagen-embedded cell, which was artificially resized by 3%, thereby simulating matrix compaction. Secondly, increasing amounts of white noise were added to the resized image in order to test the stability of both correlation methods. Relative dispersion (RD), which is defined as standard deviation divided by mean, was used as a noise level index. Finally, image resizing was followed by a horizontal translation of 60 pixels for a low and high noise example (see Table 2).

Table 2 Characteristics of validation images: CC was tested by addition of Gaussian white noise with mean 0.0 and increasing variance levels, in several cases the image was resized or translated
Full size table

3 Results

3.1 Parallel recording of matrix compaction movies

We were able to record on average 66 movies on collagen compaction by single cells in parallel at a time resolution of 15 min (n = 10 experiments). Critical issues that limited this number were the need to avoid rapid acceleration and deceleration of the microscope stage, and the time-consuming autofocus algorithm. If cells were properly focused in the initialization stage, only a small fraction (<5%) ran out of focus during 24-h compaction experiments. Under the given incubation conditions, the majority of the cells adhered firmly to the collagen matrix and demonstrated little migratory activity. Typically, cell centroids remained within a 150 μm radius of their original xy position. z position decreased slowly in many cases, reflecting compaction of the matrix in the vertical direction (data not shown).

3.2 Validation of the collagen compaction analysis

Analysis based on cross-correlation of images at a series of decomposition stages was compared with straightforward cross-correlation. This was performed on images simulating matrix compaction, subsequently followed by a challenge of increasing amounts of image noise and artificial translation.

Decomposition CC was more time-consuming than direct CC: respectively, 84 and 54 s per image pair. This was due to a larger number of cross-correlations that have to be performed, and more complex data storage and lookup operations. For low noise levels (Table 2, case c, d), both analysis methods render 100% correct displacement vectors (see Fig. 2). When noise increased, the number of vectors in the smallest decomposition stage that had insufficient cross-correlation (r < 0.5) rose. However, this effect was smaller in the decomposition CC versus direct CC. As an example, the percentage correct displacement vectors at 20.2% noise (Table 2, case f) was 77.7% in direct CC versus 89.5% in decomposition CC. Concurrently, the average displacement error at this noise level amounted to 2.1 and 0.7 pixels, respectively.

Fig. 2
figure 2

Noise sensitivity of the cross-correlation methods. Top from left to right simulated displacement field of 3% compaction (Table 2, case b), displacement field as predicted by direct CC, displacement field as predicted by decomposition CC; both cross-correlations were performed at a RD increase of 20.2% (Table 2, case f); vectors were not scaled to absolute displacements, but to optimal visual illustration; red asterisks indicate vectors based on cross-correlation values below threshold (0.5) at the highest resolution. Middle the fraction of displacement vectors in the smallest decomposition stage that is calculated correctly (error in horizontal and vertical direction ≤1 pixel). This fraction decreased at higher noise levels (relative dispersion: Table 2, case c–h). Bottom the average error between calculated and simulated displacement vectors increased for higher noise levels (Table 2, case c–h). The decomposition method resulted in a slightly lower number of incorrect vectors, and a large decrease in displacement error

Full size image

Figure 3 shows estimates for matrix compaction for these simulated deformations. The simulated compaction over a disc with 350 pixels radius was 6%. For noise levels up to 7.6% (Table 2, case d), such compaction was indeed found by both correlation methods. At higher noise levels, direct CC underestimated compaction much more severely as compared to decomposition CC. As an example, after a noise increase of 20.2%, decomposition CC estimated a compaction value of 5.9% as compared to only 4.1% for direct CC.

Fig. 3
figure 3

Calculated area at a radial distance of 350 pixels after a simulated 3% compaction. Increasing displacement errors, caused by higher noise levels, resulted in a mismatch between calculated area and simulated compacted area (94% of original). The decomposition method showed stable results at larger relative dispersion values (Table 2, case c–h)

Full size image

When an additional horizontal shift was imposed, both analysis methods correctly estimated compaction for 4.1% noise (Table 2, case i). For higher noise levels, decomposition CC succeeded in a correct area assessment. In contrast, at 20.2% noise and 60 pixels translation (Table 2, case j), direct CC predicted 1.2% rather than 6% compaction of the 350 pixels disc (see Fig. 4). This effect was not due to more incorrect displacement vectors, but could be attributed to a rise in average displacement error up to 14.6 pixels.

Fig. 4
figure 4

Sensitivity of the correlation methods to translation of the image. A horizontal displacement of 60 pixels was imposed after application of a 3% compaction and noise addition (4.1, respectively, 20.2% RD increase). Top left displacement field for high noise example (RD increase: 20.2%), as calculated by direct CC. Red asterisks indicate displacements where cross-correlation was below threshold (0.5); at these positions, a large mismatch between calculated compaction and simulated compaction occurs. Top right visual representation of matrix deformation at a radial distance of 270 pixels. Bottom calculated area remaining after a simulated compaction. At low noise levels (Table 2, case i), a horizontal shift after image resizing did not affect the calculated area; for high noise (Table 2, case j), failing correlation in the direct method results in a high offset in calculated area

Full size image

3.3 Collagen compaction by individual smooth muscle cells

Figure 5 shows an example of collagen matrix compaction by an individual smooth muscle cell. In this particular case, gross geometrical reorganization occurred in the first 6 h of the 22-h observation period. Compaction magnitude decreases with distance, and there was a time delay of around 1.5 h before compaction visible at a radius of 97 μm (110 pixels) became apparent at a distance of 308 μm (350 pixels).

Fig. 5
figure 5

Typical example of collagen matrix compaction by individual smooth muscle cell as estimated by decomposition cross-correlation, obtained at a resolution of 1 frame per 15 min. Top left displacement field represented by vectors, refined in five decomposition stages. Vectors originating from a red asterisk were created at a larger correlation scale; for the smallest window, cross-correlation value at these positions was below threshold (0.5). Top right visual representation of matrix compaction at a radial distance of 350 pixels, as estimated from a series of displacement fields, showing an initial (yellow) circle and its deformed (green) state. Bottom estimation of compacted area after 22 h at distances of 110, 230 and 350 pixels (scaling: 0.88 μm/pixel)

Full size image

4 Discussion

This study aimed at developing an automated technique for obtaining and analyzing matrix remodeling by individual cells. Emphasis was put on construction of an automatic, reliable algorithm for assessment of a detailed matrix displacement field. Especially, refinement of a cross-correlation based image analysis with a decomposition scheme was investigated. While “classic” direct CC sufficed for pairs of images with high correlation and low noise, this was no longer the case when substantial matrix remodeling occurred within the time frame between two consecutive images. This resulted in failing of CC at spots of high geometrical reorganization. However, when using decomposition CC, the gross displacements at these positions could be estimated by analysis of parent images with larger dimensions. This way, at a RD increase of 20.2% the average displacement error was lowered threefold in decomposition CC as compared to direct CC.

The method of refining displacement field accuracy with each decomposition step becomes progressively more important with larger displacements. Therefore we investigated the effect of a horizontal shift superimposed on a simulated compaction. Such shift occurred in our in vitro experiments when a group of neighboring cells pulled strongly on the matrix adjacent to a cell of interest. The shift interfered with direct CC because zero deformation was assumed when no proper correlation could be found. The result was an irregular displacement field (see Fig. 4). When using decomposition CC, on the other hand, the gross displacement was already detected in the first decomposition stage. Displacement values were then adjusted when local compaction was detected in subsequent decomposition stages. Using this strategy, a lack of local correlation at the finest resolution resulted in only a minor displacement error, with an average of 0.9 as compared to 14.6 pixels for direct CC.

Both cross-correlation methods differ not only in stability of displacement field estimation, but in efficiency of calculation time as well. Direct CC requires a large search area, i.e. the subwindows of I 2 have to extend considerably beyond the boundaries of their corresponding I 1 subwindow in order to perform a meaningful correlation. On the other hand, with decomposition CC n cc can decrease with each stage (see Table 1), since each subwindow of I 2 is repositioned according to the preliminary displacement as calculated for its parent window. Due to the increased number of calculations and data manipulations that have to be performed, analysis time was still about 50% larger for decomposition CC at the settings as stated in Table 1.

We chose a correlation threshold of 0.5 for both techniques, as well as cross-correlation windows as indicated in Table 1. These values were empirically determined as an optimum in the trade-off between calculation time, false positive displacements and overlooking local deformations. Clearly, these choices depend on the contrast and texture of the images and the nature of the deformation, and will need to be optimized for specific future experiments.

Cell traction is frequently assessed by quantification of deformations in a flexible substratum. Both 2D and 3D approaches have been used, resulting in different cellular morphologies [19, 20, 29]. While the latter seems more physiological, interpretation of the 3D experiments is more complex [35]. 2D substrates offer more straightforward tools for analyzing mechanical behavior of single cells. These materials can be enriched with fluorescent microbeads to increase image contrast. In order to achieve a high resolution in the vicinity of cells, which can change their morphology rapidly, manual tracking of specific landmarks was employed in several studies [6, 28, 30]. However, this labor-intensive method inherently limits the number of cells under study. In several studies particle tracking was performed on individual microbeads [11]. In this case, resolution depends on particle distribution, which is typically far from homogeneous in collagen substrates [23].

Several nested cross-correlation methods have been developed [35]. All of these algorithms are based on empirical determination of square size, distance for pattern search, and normalized CC threshold. In order to spatially limit the search region, a relative translational shift has been derived from image registration [23, 33]. Several schemes have been implemented to increase displacement field resolution in a step-by-step manner. In a single refinement cycle, Dembo and colleagues [26] decreased the subwindow size if significant displacement was observed. On the other hand, window size can be increased if correspondence failure occurs at high resolution [23]. The algorithm developed by Wang and coworkers [33] resembled our decomposition CC to a large extent. However, these authors aimed at construction of traction fields with a smooth nature by incorporation of filtering procedures at several stages.

Geometrical matrix reorganization provides a qualitative index of the traction forces present in the underlying material. The actual forces can be derived from a displacement field series using material properties of the ECM construct. However, for collagen scaffolds these are highly heterogenous. Frequently, local stiffness is estimated by either microneedles [5, 10, 21, 30] or optical tweezers [4, 15, 32]. Tractions can then be calculated by application of stress–strain relationships, including appropriate boundary conditions [7, 24, 25, 31, 37, 38]. However, this translation to quantitative traction forces lies beyond the scope of this article.

In conclusion, we presented integral methodology for the study of matrix remodeling. These techniques allow systematic screening of the role of matrix components such as collagen, elastin, fibronectin and laminin. Likewise, the function of stationary cells like smooth muscle cells, fibroblasts, and osteoblasts can be investigated. Our method can be applied under a wide variety of other experimental conditions. The single requirement for the image quality is a sufficiently high contrast in the material under investigation, without the neeed for laborious and potentially interfering micropatterning. Furthermore, our graphical user interface enables a flexible tuning of parameters such as number of decompositions, size of correlation window, search area and cross-correlation threshold. Using our motorized microscopy stage it is possible to patch individual images together in order to create one large field of view for the study of motile cells such as keratocytes. Finally, the method can be extended with fluorescence imaging of specific cell structures, cytokines, hormones and enzymes [15].