LSE Logo MBoC Logo

Quantifying topography-guided actin dynamics across scales using optical flow

Published Online:


The dynamic rearrangement of the actin cytoskeleton is an essential component of many mechanotransduction and cellular force generation pathways. Here we use periodic surface topographies with feature sizes comparable to those of in vivo collagen fibers to measure and compare actin dynamics for two representative cell types that have markedly different migratory modes and physiological purposes: slowly migrating epithelial MCF10A cells and polarizing, fast-migrating, neutrophil-like HL60 cells. Both cell types exhibit reproducible guidance of actin waves (esotaxis) on these topographies, enabling quantitative comparisons of actin dynamics. We adapt a computer-vision algorithm, optical flow, to measure the directions of actin waves at the submicron scale. Clustering the optical flow into regions that move in similar directions enables micron-scale measurements of actin-wave speed and direction. Although the speed and morphology of actin waves differ between MCF10A and HL60 cells, the underlying actin guidance by nanotopography is similar in both cell types at the micron and submicron scales.


Understanding the rearrangements of the cytoskeleton is essential to developing a complete picture of the dynamic forces involved in cellular processes such as migration, division, and differentiation. Cytoskeletal dynamics, and in particular actin dynamics, have been shown to be important for the growth of cell junctions and focal adhesions (Maruthamuthu et al., 2010) and for immune-cell activation (Lam Hui et al., 2014). The formation of actin waves through directional polymerization and depolymerization of filaments drives many types of cell migration (Müller and Sixt, 2017) and has been associated with the establishment of polarity in a variety of cell types (Inagaki and Katsuno, 2017).

Forces from the extracellular environment are an important modulator of actin dynamics. Physical and chemical characteristics of the extracellular environment, such as rigidity, biochemical composition, and topography, have been shown to influence actin dynamics and associated cell behavior (Discher et al., 2005; Doyle et al., 2009; Lu et al., 2012; Petrie and Yamada, 2015). One mechanism for this modulation is mechanosensing via focal adhesions (Ketchum et al., 2018). In addition, actin waves respond when cells encounter obstacles (Weiner et al., 2007). It has been established that ridges of width comparable to fibers in the extracellular matrix (ECM) can alter actin dynamics significantly (Driscoll et al., 2014; Sun et al., 2015; Ketchum et al., 2018) and bias the localization of focal adhesions (Ventre et al., 2014; Sun et al., 2018). Thus, in vivo, the topography of the ECM, such as collagen networks (Wolf and Friedl, 2011; Charras and Sahai, 2014), is likely to modulate actin dynamics.

Periodic nanotopographic surfaces provide the opportunity to obtain systematic data on the modulation of such intracellular dynamics. In prior work, we have shown that actin waves can be nucleated near, and guided along, periodic nanotopography, in a phenomenon termed esotaxis. Actin-wave guidance has been observed in cell types that exhibit distinct physiological functions and migration phenotypes, including Dictyostelium discoideum (Driscoll et al., 2014; Sun et al., 2015), neutrophil-like HL60 cells (Sun et al., 2015), B cells (Ketchum et al., 2018), and breast cancer cell lines (Chen et al., 2019). However, there are clear differences in the responses of each of these cell types to nanotopography. For example, although both D. discoideum and HL60 cells exhibit esotaxis, these two types of cells have been found to move preferentially in different directions on specific nanoscale asymmetric sawtooth textures (Sun et al., 2015). Furthermore, different breast cancer cell lines preferentially move in different directions on asymmetric sawtooth nanotopography (Chen et al., 2019).

Here we introduce a method for performing quantitative measurements of the influence of nanotopography on intracellular dynamics at both the submicron and the micron scales. This approach enables the detection of subtle differences in cytoskeletal dynamics and allows for in-depth analysis of both the differences and the similarities of these dynamics across cell types and phyla. Our method of quantification of actin dynamics across scales is based on optical flow, an image-analysis technique developed in the fields of robotics and navigation control that uses changes in pixel intensities to detect motion in image sequences (Horn and Schunck, 1981; Lucas and Kanade, 1981). Because of the popularity of particle image velocime­try (PIV), optical flow has seen limited use on biological images. However, PIV is poorly suited for the variety of features that can be exhibited in fluorescence images of amorphous concentration fields. Indeed, a recent study indicated that optical flow may be better suited for analysis of fluorescence images, as it provides a more accurate estimate of ground-truth flow fields (Vig et al., 2016). Here, we use optical flow to measure the dynamics of actin polymerization with submicron precision, and we further expand the utility of optical flow by introducing modeling and fitting approaches to the analysis of optical-flow vector fields. Clustering of the optical-flow data further allows us to quantify actin dynamics on the micron scale. Thus, this optical–flow-based analysis enables the identification of similarities and differences between esotaxis in neutrophillike HL60 cells and human breast epithelial MCF10A cells across length scales.


Esotaxis has been observed in a wide range of cell types that are known to respond to their in vivo microenvironment through processes such as directed migration or immune-system activation (Driscoll et al., 2014; Sun et al., 2015; Ketchum et al., 2018). More recently, esotaxis has also been observed in epithelial cells, which are less motile (Chen et al., 2019). Here we contrast the actin dynamics of epithelial MCF10A cells with those of neutrophillike HL60 cells.

LifeAct–GFP-labeled epithelial MCF10A cells were plated on a 900 µm × 900 µm region patterned with parallel nanoridges with a spacing of 1.5 µm, as well as on the surrounding flat region. Confocal imaging near the surface revealed distinct actin morphologies on nanoridges as compared with the flat region (Figure 1). On the flat region the phenotype is a common one for these cells on such surfaces, with a broad lamellipodium at the cell front and stress fibers throughout the cell body (Figure 1A). In contrast, MCF10A cells on nanoridges exhibit actin streaks aligned with the ridges throughout the cell area (Figure 1B). The local nature of the response of actin to surface texture is illustrated in Figure 1C, which shows a cell that is partially on the nanoridges and partially on the flat region. On the nanoridged region, the cell shows the same actin streaks as a cell that lies fully on a ridged surface, whereas the same cell maintains a broad lamellipodium on the flat region.


FIGURE 1: Surface topography prompts distinct actin morphology. The actin cytoskeleton of an epithelial MCF10A cell on (A) a flat surface has an actin morphology that is distinct from that of a cell on (B) a nanoridged surface. (C) A cell that is partially on a nanoridged region and partially on a flat region exhibits local actin morphologies that are driven by the underlying topography. The blue line in C indicates the boundary between the flat region and the nanoridged region. All scale bars are 10 µm. See also Supplemental Movie S1.

Movie S1: Time‐lapse images of the MCF10A LifeAct‐EGFP cells shown in Figure 1. Images were taken every 10 s and are shown at 20 frames/s; the clock in the upper left corner indicates min:s. The scale bar indicates 10 μm.

Ridged and flat regions also engender distinct actin dynamics. Kymographs can be used to visualize dynamics in a region of interest in one spatial direction over time. The left side of Figure 2A compares an MCF10A cell on a flat region with one on a nanoridged region. The cell on the nanoridged region shows actin streaks that are characteristic of esotaxis. Actin kymographs from two perpendicular regions (Figure 2A) in an MCF10A cell on a flat surface show oscillatory dynamics in all directions at the cell boundary. These oscillations in the kymographs indicate the presence of fanlike protrusions and retractions across each region over 30 min (Figure 2B). In contrast, on the nanoridged region, the actin dynamics parallel and perpendicular to the ridges are different (Figure 2C). Parallel to the nanoridges, MCF10A cells show oscillatory actin dynamics. As shown in the bottom left of Figure 2C, a representative kymograph of a region perpendicular to the ridges shows actin structures that persist for tens of minutes and do not move perpendicular to the ridges. This behavior is typical for kymographs perpendicular to the ridges, although perpendicular motion is observed in some cases as discussed below.


FIGURE 2: Surface topography leads to distinct actin morphodynamics. (A) Optical micrographs of MCF10A cells in (left) bright-field and (right) fluorescence on (top) a flat region and (bottom) a nanoridged region. All scale bars are 10 µm. Kymographs for the areas denoted in A are shown in B for the flat region and in C for the nanoridges. D–F are the same as A–C, respectively, but for HL60 cells. The dynamics in A and D are shown in Supplemental Movies S2 and S3, respectively.

Movie S2: Time‐lapse images of the MCF10A LifeAct‐EGFP cells shown in Figure 2. Images on the left are bright field images corresponding to the actin shown on the right.Images were taken every 10 s and are shown at 20 frames/s; the clock in the upper left corner indicates min:s. The scale bar indicates 10 μm.
Movie S3: Time‐lapse images of the HL60 Actin‐YFP cells shown in Figure 2. Images on the left are bright field images corresponding to the actin shown on the right. Images were taken every 2 s and are shown at 20 frames/s; the clock in the upper left corner indicates min:s.The scale bar indicates 10 μm.

The behavior of motile, neutrophil-like HL60 cells on flat and nanoridged regions is illustrated in Figure 2, D–F. Figure 2D illustrates the regions from which kymographs were generated. In HL60 cells on flat regions, actin is concentrated near the cell front. This localization is largely preserved on the ridged surfaces, although the morphology of the actin changes such that streaks of actin are aligned with the ridges. On flat surfaces, the HL60 cells show regions of protrusions and retractions (Figure 2E) similar to the actin dynamics seen in MCF10A cells (Figure 2B). We note that protrusions occur on the scale of seconds in HL60 cells and on a scale of minutes in MCF10A cells. Kymographs of the HL60 cells in the direction parallel to the nanoridges show protrusive dynamics, although often in the form of a single persistent wave (Figure 2F, top), in contrast to the oscillatory behavior seen on flat surfaces (Figure 2E). A representative kymograph of an HL60 cell in the direction perpendicular to the nanoridges shows streaks (Figure 2F, bottom) that indicate that actin waves do not move perpendicular to the ridges, but the streaks are shorter in duration than those in a typical MCF10A cell (Figure 2C). This behavior is typical for kymographs of actin in HL60 perpendicular to the ridges. Unlike in MCF10A cells, in which actin streaks on the ridges localize throughout the cell (Figure 2C, bottom), in the HL60 cells the streaks occur near the cell front (Figure 2F, bottom). Groups of actin streaks propagate together at the front of the HL60 cells, suggesting that there may be large-scale organization of actin dynamics (spanning many ridges) in these cells. It is unclear whether there is large-scale organization of actin dynamics in the MCF10A cells.

As shown in Supplemental Movies S2 and S3, the full range of actin dynamics is more complex than is revealed by kymographs. MCF10A cells on the ridged regions exhibit actin dynamics throughout the substrate contact area, whereas actin dynamics on flat surfaces are largely confined to the cell boundary. These movies show that, in both cell types, nanoridges stimulate reproducible, dynamic, linear actin structures.

Time-lapse fluorescence images of actin waves are difficult to interpret by visual inspection or kymographs alone, because the observed dynamics arise from a complex spatio-temporal concentration field. To measure these wavelike dynamics quantitatively, we must first define a wave (size and shape) and then capture its propagation (splitting, recombination, and changes in direction). Here, we address these challenges by introducing an automated approach to quantify actin-wave dynamics across length scales for unbiased comparison in different cell types and extracellular environments.

Our method is based on a computer-vision algorithm from robotics and navigation control called optical flow (Horn and Schunck, 1981; Lucas and Kanade, 1981), which provides pixel-based information about the direction and magnitude of intensity flux in a series of time-lapse images. Fields of optical-flow vectors are calculated by integrating changes of intensity in space and time, as shown schematically in Figure 3. For example, two images of a migrating HL60 cell taken 8 s apart are shown with changes in time highlighted by a green-to-magenta montage (Figure 3A). The magenta region indicates growth of the actin front (which, as expected, occurs at the leading edge of the cell), and the green region indicates a decrease in actin intensity.


FIGURE 3: Optical-flow calculations capture the dynamics in movies of actin fluorescence. (A) Two frames of a representative HL60 cell obtained 8 s apart and a merged image show the dynamics of the cell’s behavior over time. The schematic in B illustrates how the procedure used to carry out optical-flow calculations combines the spatial gradient of an image (left) and the difference image/temporal gradient (right) to yield the optical-flow vector field (center). These calculations are applied to the images in A and shown in the images of C. The spatial gradient field (left) and temporal gradient (right) result in the output optical-flow vector field (center). Blue pixels in the right panel of C indicate a positive change (increase) in the pixel brightness from the first frame to the second frame, and red pixels in the right panel indicate a negative change (decrease) in the pixel brightness from the first frame to the second frame. All scale bars are 5 µm.

The general objective of calculating optical flow is to solve for the unknowns Δx and Δy in


where I(x,y,t) represents the actin fluorescence intensity at frame t. The intensity I that exists at point (x,y) at time t translates to a new point (x + Δx,yy) at some future time t + Δt. Expanding about small Δx and Δy and neglecting second-order derivatives yields the master optical-flow equation


This governing equation is underdetermined, and so the Lucas–Kanade optical-flow constraint (Lucas and Kanade, 1981) was applied to calculate flow fields. This constraint prescribes that all pixels in a small window centered at (x,y) each have the same translational optical-flow vector. The equation can then be solved using the least-squares criterion (an explicit derivation is given in the Materials and Methods) to yield the intensity flow, (Figure 3B, center panel). Solving for requires use of the negative spatial gradient, (Figure 3B, left panel), which forms a vector field oriented away from regions of highest local intensity, and the time derivative, ∂I/∂t (Figure 3B, right panel) as shown in Eq. 2. Using the pair of images from Figure 3A as a representative example, the spatial and temporal gradients are used to calculate the optical-flow vector field, which approximates the flow of actin between the two frames. In this example, the vector field captures the translational motion on the leading edge of the cell (Figure 3C).

Optical-flow measurements of actin intensity translation enable the quantification of the pixel-scale response of actin dynamics to nanoridge topographies (Figure 4). The green-to-magenta montages of representative HL60 and MCF10A cells show dynamic and protrusive actin behavior at the leading edge of the cell (Figure 4A). Coloring the calculated flow fields based on direction relative to the nanoridges (Figure 4B) reveals the clear bias of actin wave guidance in the direction parallel to the ridges, which is consistent with the qualitative features of the montage images in Figure 4A. Measurements of the optical-flow directions on all HL60 and MCF10A cells on both flat and ridged surface topographies are shown in the histograms of Figure 4C. In both cell types, the cumulative distribution of flow in cells on flat surfaces shows no appreciable bias in any direction. However, cells on ridges exhibit a clear preference for flow along the ridge direction.


FIGURE 4: Using optical flow to measure pixel-scale guidance. (A) Representative merged time-lapse images of an MCF10A cell (top, 100 s apart) and an HL60 cell (bottom, 8 s apart). (B) Optical-flow vector fields colored by direction relative to the horizontal ridges. Blue indicates motion aligned with the ridges and yellow indicates motion perpendicular to the ridges. All scale bars are 5 µm. The dynamics in B are shown in Supplemental Movies S4 and S5. (C) Cumulative distributions of the directions of optical-flow vectors for multiple HL60 and MCF10A cells on flat and ridged surfaces; all cells are weighted equally in the distribution. (N = 4 MCF10A cells on flat surfaces from three independent experiments, N = 19 MCF10A cells on ridges from four independent experiments, N = 14 HL60 cells on flat surfaces from two independent experiments, and N = 17 HL60 cells on ridges from three independent experiments.) (D) The distribution of angles can be fit to a mixture of two von Mises distributions with five fitting parameters: θµ (primary direction of motion), κ (inversely related to distribution width), and three coefficients indicating the component of motion in the θµ direction, the component in the θµ + 180° direction, and the component that is uniform. (E) In both MCF10A and HL60 cells, the distribution width, parameterized by κ, shows significant differences (p = 0.0001583 and p = 0.0040) on flat versus ridged surfaces. (F) For each cell, the mean direction of motion (angular axis) is plotted vs. κ (radial axis). Values of κ less than 2 (indicated by the dashed line) indicate cells with direction distributions that are statistically indistinguishable from a uniform distribution. We note κ = 19.6 for one HL60 cell, and this point is not visible in this figure.

Movie S4: Time‐lapse images of the optical flow results for the example MCF10A cell shown in Figure 4B. Actin fluorescence images are overlaid with color‐coded optical‐flow vectors. Blue vectors are pointing left and right, and yellow vectors are pointing up and down.Images were taken every 10 s and are shown at 12 frames/s; the clock in the upper left corner indicates min:s.
Movie S5: Time‐lapse images of the optical‐flow results for the example HL60 cell shown in Figure 4B. Actin fluorescence images are overlaid with color‐coded optical‐flow vectors. Blue vectors are pointing left and right, and yellow vectors are pointing up and down.Images were taken every 10 s and are shown at 12 frames/s; the clock in the upper left corner indicates min:s.

We note that the images for the MCF10A and HL60 were obtained with different objectives (100× and 60×, respectively) and at different acquisition rates (10 and 2 s between frames, respectively). Nevertheless, the optical-flow algorithm performs well on each set of data, emphasizing the broad applicability of this approach. The two cell types also use different actin-labeling approaches (LifeAct for the MCF10A and Actin-YFP for the HL60), which may have different overexpression effects on actin dynamics (Melak et al., 2017), but in both cases produce reliable optical-flow results that are indicative of nanotopography-guided actin dynamics.

For further quantification, we fit the distribution of flow directions from each cell to a bimodal von Mises model with a constant offset (Materials and Methods). The distribution used consists of a uniform component and two peaked components that are 180° apart. The five parameters of the bimodal model are illustrated in Figure 4D. The angle θµ indicates the direction of the main component, and 1/κ is proportional to the width of the distribution. The values of κ on ridged regions are significantly higher than those than on flat regions for both the MCF10A and the HL60 cells (Figure 4E, p = 0.0001583 and p = 0.0040), indicating that the ridges strongly guide the actin flows in a bidirectional manner. A comparison of κ and θµ shows that cells with a bidirectional actin flow (i.e., high κ values) are more likely to be guided along the ridge direction (Figure 4F).

Although the optical-flow vector field indicates preferred directions of actin flow, it does not yield propagation speeds of actin polymerization waves directly. The magnitude of an optical-flow vector incorporates both the shift of actin in space and the change in actin intensity over time. This submicron-scale (i.e., pixel-scale) flux of intensity does not translate directly into characteristics of the dynamics that are notable on the micron-scale (i.e., tens of pixels), such as the organization of waves across ridges seen in neutrophils in Figure 2F or the speed of wave propagation. To quantify these micron-scale characteristics, we combined similarly oriented optical-flow vectors into clusters (Figure 5A), which were then tracked over time. To ensure that we tracked robust clusters, we applied additional constraints, such as requiring sufficiently large intensity changes (see Materials and Methods). The result of this clustering was the identification of broad regions of actin that moved collectively (Figure 5B).


FIGURE 5: Clustering of optical-flow vectors to measure micron-scale dynamics. (A) Similar flow vectors are grouped into clusters. (B) Clusters contain optical-flow vectors with a wide array of orientations, resulting in micron-scale structures. All scale bars are 5 µm. (C, D) Particle-tracking algorithms are applied to the tracked clusters. The cluster tracks are consistent with the motion at the leading edge of the actin waves seen in the kymographs in both HL60 (C, left) and MCF10A (D, left). Panels to the right of each kymograph show the same track in the 2D context of the cells; clusters found throughout the cell over time are indicated by colored regions. The dynamics in C and D are shown in Supplemental Movies S6 and S7, respectively. (E) The cluster tracks are used to determine speed distributions of actin waves on the ridges in the MCF10A and HL60 cells, which show no significant difference (p = 0.5529 and p = 0.0586) between flat surfaces and ridges. (F) Cumulative angle distribution of cluster track directions for multiple HL60 and MCF10A cells on flat and ridged surfaces; all cells are weighted equally in the distribution. (N = 4 MCF10A cells on flat surfaces from three independent experiments, N = 17 MCF10A cells on ridges from four independent experiments, N = 9 HL60 cells on flat surfaces from two independent experiments, and N = 16 HL60 cells on ridges from three independent experiments.)

Movie S6: Time‐lapse images of the HL60actin track shown in Figure 5C. Actin images are overlaid with clusters regions and cluster tracks. The clusters are color‐coded to match the times indicated in Figure 5. Images were taken every 10 s and are shown at 12 frames/s; the clock in the upper left corner indicates min:s.The scale bar indicates 5 μm.
Movie S7: Time‐lapse images of the MCF10A actin track shown in Figure 5D. Actin images are overlaid with clusters regions and cluster tracks. The clusters are color‐coded to match the times indicated in Figure 5. Images were taken every 10 s and are shown at 12 frames/s; the clock in the upper left corner indicates min:s.The scale bar indicates 5 μm.

We applied peak-finding and tracking algorithms (Blair and Dufresne, 2018) to follow the locations of maximum alignment of these optical-flow clusters on the micron scale. Although the optical-flow results shown in Figure 4 follow motion on the pixel scale, by following the peak alignment of the optical flow, we are able to track larger coordinated clusters of actin. Our tracking is consistent with the actin dynamics seen in kymographs, as shown by representative kymographs of MCF10A (Figure 5C) and HL60 (Figure 5D) cells overlaid with the tracked location of actin waves. Unlike kymographs, which are sensitive to motion along a chosen direction, tracks of clustered flow vectors reveal the micron-scale motion of actin in two dimensions. The benefits of our approach are illustrated in Figure 5C for an MCF10A cell. For the initial frames, the kymograph indicates a stationary actin structure (Figure 5C, left), but when the actin dynamics are viewed in two dimensions (Figure 5C, right), it is evident that in the early frames this wave is moving perpendicular to the ridges, a motion that cannot be captured in this one-dimensional kymograph. Thus, the combination of optical flow, clustering, and tracking allows us to follow actin waves without being limited to tracking only motion that occurs along a straight line.

The speeds of the tracked actin clusters (Figure 5E) are similar to speeds derived from actin kymographs (Figure 5, C and D, and Supplemental Figure S1), despite an approximately order-of-magnitude difference in speed between the two cell types that is consistent with their distinct in vivo functions and with previously reported cell-migration speeds (Meyer and Howard, 1987; Lee et al., 2016). For both cell types we find no significant difference (p = 0.5529 and p = 0.0586) between actin-wave speeds on flat and ridged surfaces (Figure 5E), implying that topography steers actin dynamics but does not alter wave speeds. On flat surfaces, the directions of the clusters are distributed uniformly for MCF10A cells but show distinct peaks in multiple directions for HL60 cells (Figure 5F). This observation is consistent with the polarized character of actin in several of the HL60 cells on flat surfaces, corresponding to κ values greater than 2 in Figure 4, E and F.


Extracellular texture, which is an important component of the 3D, in vivo environment, is capable of spatially patterning actin and modulating actin dynamics. Using nanoridge structures in conjunction with optical-flow approaches, we are able to probe and quantify this intracellular response to extracellular textures in a systematic manner.

Previous studies of D. discoideum (Driscoll et al., 2014; Sun et al., 2015), B cells (Ketchum et al., 2018), and tumor-associated fibroblasts (Azatov et al., 2017) showed similarity in actin response to texture, which suggests that guidance of actin driven by texture (esotaxis) is broadly conserved across cell types. Controlled textures are thus a useful model microenvironment for the systematic, reproducible, and quantitative study of intracellular dynamics and force regulation. Here we demonstrated the analysis of time-lapse images of two cell types that have distinct physiological function. Neutrophil-like HL60 cells are polarized and highly motile, and respond to a variety of cues as they fulfill their role in the immune system. Epithelial MCF10A cells, however, have a nonmotile physiological function. Nevertheless, both cell types show clear, and quantitatively similar, actin dynamics in response to surface textures. Consistent with our prior results (Driscoll et al., 2014; Azatov et al., 2017), we find that nanoridges lead to persistent streaks of actin that are not seen on flat surfaces.

Optical flow enables the quantification of both the reproducible streaks of actin seen on nanoridged surfaces and the more chaotic actin waves seen on flat surfaces. The latter waves are typically much wider than guided actin waves. On flat surfaces the waves often change direction and can also grow wider and split. Such motion phenotypes are not easily captured with standard techniques such as kymographs. Optical flow enables us to follow these dynamics and thus yields insights beyond those derived from kymograph-based techniques.

We note that optical flow is suitable for comparisons of systems imaged under different conditions (e.g., 60× vs. 100× objectives), enabling comparisons of widely varying cell sizes and migration speeds. The use of varying acquisition rates (i.e., 2 and 10 s) in this work was based on the differences in cell migration speeds of the two cell lines studied. In general, optical flow requires a frame rate such that changes in fluorescence intensity between frames are small but larger than noise. Our use of the Lucas–Kanade optical-flow constraint also makes the assumption that there is a smoothness to the flow field over a certain neighborhood, which is a length-scale parameter in the optical-flow analysis. This assumption is met by a wide variety of biological imaging data sets, and thus the use of our optical-flow approach is not limited to actin dynamics. Optical flow could provide insights into the motion of other cytoskeletal proteins, such as tubulin, or into the dynamics of other fluorescent markers that exhibit a spatially and temporally changing intensity field. Our use of clustering to study larger-scale actin dynamics could similarly be adapted to other fluorescent markers under the assumption that there are larger-scale dynamics that move together in similar directions. In this work, we used a spinning-disk confocal microscope for image acquisition, but our analysis pipeline would also be appropriate for other imaging techniques, such as epifluorescence. When working with other imaging modalities or fluorescent markers, the size of the Lucas–Kanade neighborhood and the threshold for vector reproducibility (see Materials and Methods and Supplemental Figure S2) can be adapted to only include robust results in further analysis.

Using submicron-scale optical flow and associated micron-scale analysis, we have shown that both MCF10A and HL60 cells have actin flows that are biased along nanoridges. By clustering similarly oriented optical-flow vectors, we are able to measure the speed of actin waves within the cell. The measured speeds are comparable to speeds calculated from actin kymographs. Optical-flow analysis allows us to determine that the speeds do not differ significantly on flat versus ridged regions. This finding indicates that nanotopography guides, but does not fundamentally alter, the speed of actin dynamics. We measure actin-wave speeds on the order of 1 µm/min in the MCF10A cells, consistent with previously reported cell migration speeds of approximately 0.5 µm/min (Lee et al., 2016). In the HL60 cells we find actin speeds ranging from approximately 8 to 14 µm/min, consistent with the 8 µm/min speed for cell migration previously reported (Meyer and Howard, 1987).

Fitting the optical-flow vectors to a bimodal von Mises distribution enables quantification of the differences in the directionality of actin flows on flat and ridged surfaces in both cell lines. The fit parameters also show differences in actin polarization in these two cell lines. HL60 cells occasionally exhibit coordinated and directed actin flows even on flat surfaces, whereas MCF10A cells on flat surfaces show uniform direction distributions of actin waves. On the micron scale, actin-wave tracks from individual HL60 cells on flat surfaces generally polarize and have a preferred direction, consistent with the behavior of immune cells, which tend to polarize and migrate in a directed manner. Tracks from MCF10A cells on a flat surface, on the other hand, are more directionally uniform for each cell. In both cell types, ridges guide actin waves in a bidirectional manner.

The quantitative actin responses in MCF10A and HL60 cells support a model in which surface texture provides a symmetry-breaking cue that leads to nucleation of actin polymerization. Under flat tissue-culture conditions, which lack symmetry-breaking cues, actin polymerization relies on spontaneous nucleation or edge effects (Begemann et al., 2019). Edge effects may lead to morphological features such as the lamellipodia seen in HL60 cells on flat surfaces in Figure 2B. By changing the landscape on which nucleation occurs, surface texture can lead to actin polymerization in other locations of the cell, such as the persistent streaks seen in the center of MCF10A cells on ridges in Figure 2A.

There are multiple mechanisms by which cells may respond to local forces and geometry (Vogel and Sheetz, 2006), including sensing mechanisms that can respond to membrane curvature on a variety of scales (McMahon and Boucrot, 2015). Our finding that nanoridges change the direction, rather than the speed, of actin waves suggests that growth of existing actin filaments away from the surface is the rate-limiting step in actin polymerization wave propagation. In some cases, sensing mechanisms may rely on the preferential formation of focal adhesions. This hypothesis is consistent with previous results on focal adhesion localization and orientation in response to surface texture (Ventre et al., 2014; Sun et al., 2018). Although MCF10A cells form strong focal adhesions that may align with texture cues (Sun et al., 2018), the HL60 cells form weaker adhesions, and the previously studied D. discoideum cells (Driscoll et al., 2014; Sun et al., 2015) are not known to form integrin-mediated focal adhesions. Thus, the dominant mechanism of surface texture response likely depends on both the cell type and the extracellular environment.

Known surface-sensing mechanisms also include cytoskeletal components such as septins, which respond to micron-scale curvatures (Bridges et al., 2016), and BAR domains, which sense nanoscale curvature (Zimmerberg and McLaughlin, 2004). Proteins with BAR domains have been linked to actin assembly (Graziano and Weiner, 2014) as well as to key components of actin-regulating pathways, such as WAVE and Rac (Miki et al., 2000; Habermann, 2004). Recent work has suggested that nucleation of new actin filaments is enhanced by nanotopography. Specifically, curved nanopillars activate the nucleation-promoting factors Arp2/3 and N-WASP through enhanced binding of an F-BAR domain containing protein (Lou et al., 2019). Additionally, evidence suggests that topography is capable of shifting multiple gene-expression pathways (Dalby et al., 2005), which implies that longer-term exposure to topography may mediate additional surface-sensing pathways. As in vivo microenvironments contain a variety of textures, it is likely that multiple mechanisms respond to distinct features of extracellular texture, and future work on the response of actin regulators to controlled topographies such as those investigated here will help elucidate the contributions of distinct signaling pathways in topography-guided actin dynamics.

Although the systematic modulation and interrogation of all possible molecular factors of esotaxis is beyond the scope of this article, our analysis yields two remarkable constraints on the molecular sources of esotaxis. First, the speed of actin waves is not altered by esotaxis. Second, the directional guidance provided by nanotopography is comparable in the two cell types investigated, despite their disparate functions and migratory phenotypes. Quantitative analysis of esotaxis as a physical phenotype could yield crucial prognostic disease insights, especially in the case of cancer, in which changes in the texture of the microenvironment correlate with disease progression.


Cell culture and imaging

HL60 YFP-actin cells were a gift from Orion Weiner of the University of California, San Francisco. The cells were cultured in RPMI 1640 medium, Glutamax (Life Technologies) supplemented with 10% heat-inactivated fetal bovine serum (Gemini Bio). Cells were passaged every 2–3 d and kept between 3 × 105 and 1 × 106 cells/ml. For differentiation, cell media was additionally supplemented with 1.3% dimethyl sulfoxide Hybri-Max (Sigma Aldrich) for 5 d before imaging. Actin dynamics of HL-60 YFP-actin cells were observed by confocal fluorescence and bright-field time-lapse imaging using a PerkinElmer spinning-disk confocal microscope with a water immersion 60× objective (0.21 µm/pixel). Images were recorded every 2 s. We note that this method of plating resulted in the imaging of some multicellular clusters of HL60 cells; these clusters were removed from further analysis.

Preparation for imaging included a 10 µg/ml coating of fibronectin on the substrates. Cells were plated and allowed to settle. After approximately 30 min, N-Formyl-Met-Leu-Phe (fMLF; Sigma Aldrich) was added to 1 μM. HL-60 cells were imaged on flat resin and ridged nanotopographies beginning between 10 and 15 s after fMLF stimulation. All images analyzed in this work were obtained after fMLF stimulation.

MCF10A LifeAct-EGFP cells were a gift from Carole A. Parent (National Cancer Institute, Bethesda, MD). These cells were cultured in DMEM/F12 media supplemented with 5% horse serum, 10 μg/ml insulin (ThermoFisher Scientific), 10 ng/ml EGF (Peprotech, Rocky Hill, NJ), 0.5 μg/ml hydrocortisone, and 100 ng/ml cholera toxin (both Sigma, St. Louis, MO). The media were additionally supplemented with 2 μg/ml puromycin dihydrochloride (ThermoFisher Scientific) to select for EGFP-positive cells. Before imaging, cells were plated on a nanoridged surface coated with collagen IV and were allowed to adhere to the surface overnight. Actin dynamics were studied by confocal fluorescence and bright-field, time-lapse imaging using a PerkinElmer spinning-disk confocal microscope with a 100× objective (0.14 µm/pixel). Images were recorded every 10 s.

For both cell types, data were collected using PerkinElmer’s Volocity software (version 6.4.0). The spinning-disk confocal microscope was equipped with a Hamamatsu ImagEM X2 EM-CCD camera (C9100-23B), which recorded 12-bit images. Cells used in this study tested mycoplasma negative using the MycoAlert (Lonza) testing system.

Surface fabrication

The nanotopographies were designed and fabricated using multiphoton absorption polymerization (MAP), as described previously (Sun et al., 2018). A drop of prepolymer resin (1:1 wt/wt Tris [2-hydroxy ethyl] isocyanurate triacrylate [SR368]: ethoxylated (6) trimethylolpropane triacrylate [SR499], both from Sartomer; 3% Lucirin TPO-L [BASF]) was sandwiched between a coverslip and a plasma-treated microscope slide that had been functionalized with acrylate groups (LaFratta et al., 2007; Li and Fourkas, 2007; Sun et al., 2018). The coverslip was mounted onto the stage of an inverted microscope (Zeiss Axiovert 135). A beam of 150-fs pulses centered at 800 nm from a Ti:sapphire oscillator (Coherent Mira 900) was directed into the microscope and through a high-numerical-aperture objective (Zeiss alpha-Plan Fluar 100×; NA 1.45). The stage motion and shutter were controlled using a program written in LabVIEW (National Instruments). Once the pattern was fabricated, the sample was developed in ethanol and baked at 110°C for 1 h.

A replica molding approach was then used to mold the chemically functionalized pattern (Sun et al., 2018). This step included the initial casting of a hard-poly(dimethylsiloxane) (h-PDMS) layer (1000 rpm, 40 s) to better resolve nanoscale features of the topographical pattern. This layer was allowed to sit on the pattern at room temperature for 2 h and was then baked at 60°C for an additional 1 h. Finally, Sylgard 184 was mixed (10:1 elastomer base:curing agent) and poured onto the initial h-PDMS layer. The PDMS was cured at 60°C for 1 h 10 min.

The MAP-fabricated structure was then reproduced through a soft-lithographic technique. A drop of the aforementioned resin was sandwiched between the PDMS mold and an acrylate-functionalized coverslip and was then exposed to UV radiation from a lamp for a desired amount of time. After the resin cured, the coverslip was peeled off the mold. This process was repeated many times to produce enough replicas to perform the necessary experiments. The replicas were soaked in ethanol for at least 4 h and subsequently baked/dried in an oven at 110°C for 1 h. Samples used to study MCF10A cells were also soaked in UltraPure water (ThermoFisher) for approximately 12 h.


Kymographs were created in MATLAB by manually selecting a rectangular region in an actin image. Fluorescence intensities inside the region were averaged across the short axis of the region; this process was repeated for each image in the time-lapse sequence, and the resulting intensity data were combined into the kymographs shown in Figure 2 and Supplemental Figure S1.

Optical flow

The Lucas–Kanade optical-flow method (Lucas and Kanade, 1981) was used to capture the direction and strength of intensity flow of fluorescent actin and to produce vector fields indicating actin motion. The optical-flow vector field of an image series is the field of apparent translation in the image plane, as is shown schematically in Figure 3B. Calculating the optical flow for two adjacent two-dimensional images in an image series requires solving for the unknowns Δx and Δy in Eq 1, as described above. The Lucas–Kanade method uses a least-squares regression approach to solve for the best optical-flow vector on a pixel-by-pixel basis under the assumption that all pixels within a “neighborhood” move in a similar direction (Lucas and Kanade, 1981). If solving for the optical-flow vector of some point p with coordinate (xp,yp), the master optical-flow equation requires that the optical flow vector of point p and all points in the neighborhood about p (points 1,2,…,p,…, N) follow the underdetermined relation


Under the Lucas–Kanade assumption, the vector for point p is assigned to all points in the neighborhood


The least-squares solution to equations of this form, , is

. Furthermore, we implement a scheme using a Gaussian weight matrix centered at point p to ensure that pixels near p have more influence over the result. The equation then becomes


where w is a centered Gaussian matrix. The equation then takes the form and has solution


Optical-flow reliability is defined as the smallest eigenvalue of the ATwA matrix (Simoncelli et al., 1991; Barron et al., 1994) and was used to remove flow vectors that were created by noise or an ill-defined least-squares calculation. The threshold used can be adapted to best suit the experimental data and scientific questions of interest by only keeping the most reliable vectors while measuring the motion of more regions of the cell. The optical-flow weight matrix for MCF10A cells was a 27 × 27 pixel Gaussian with a SD of 4.5 pixels (0.63 µm). The optical-flow weight matrix for HL60 cells was 19 × 19 pixel Gaussian with a SD of 3 pixels (0.63 µm). With other imaging modalities, magnifications, or fluorescent markers, it may be appropriate to change the size of this weight matrix based on the size of features of interest and noise in the image.

von Mises model of flow distribution

Optical-flow distributions were modeled with a modified bimodal von Mises distribution (von Mises distributions are a continuous and differentiable analogue of normal distributions on a circle with similar statistical properties). The model was defined as


where VM (θ | θμ,k) is the von Mises distribution


and is the modified Bessel function of the first kind. The maximum likelihood estimates of the parameter k were used for statistical analyses.

Cluster-tracking analysis

Regions of actin fluorescence were clustered using the direction of optical-flow vectors together with an optical-flow reliability threshold and by requiring that actin intensity change over time (see Supplemental Figure S2 for a visualization of this workflow). The dot products between optical-flow vectors around a point p (i.e., vectors v1,v2,…,vp-1, vp+1,…,vN) were calculated and accumulated using a Gaussian weighting scheme to a single scalar alignment metric. The alignment metric is defined as


where w is a renormalized centered Gaussian matrix with a center manually set to 0. This calculation was carried out for each pixel.

To require that the actin intensity change over time, a mask of the thresholded difference image between subsequent frames was calculated. For every pair of adjacent frames, It and Itt, the resulting mask took value I where Itt > It and 0 otherwise. For our analysis, Δt = 30 s for MCF10A and 6 s for HL60.

To calculate the final clustered regions, the alignment metric ap, optical-flow reliability λp, and difference-image mask were multiplied in an element-wise manner to create a final cluster image. The cluster image was inputted into a peak-finding algorithm to locate peaks in the resulting intensity profile, and the Crocker-Grier particle-tracking algorithm (Crocker and Grier, 1996; Blair and Dufresne, 2018) was used to track coordinates of the resulting peaks over time.

The clustering weight matrix for MCF10A was a 27 × 27 pixel Gaussian with a SD of 4.5 pixels (0.63 µm). The clustering weight matrix for HL60 cells was 19 × 19 pixel Gaussian with a SD of 3 pixels (0.63 µm). The diameter of the peaks used in pkfnd.m (Blair and Dufresne, 2018) was 15 pixels (2.1 µm) for MCF10A cells and 10 pixels (2.1 µm) for HL60 cells. The maximum displacement used in track.m (Blair and Dufresne, 2018) was 11.5 pixels (1.54 µm) for MCF10A cells and 7 pixels (1.47 µm) for HL60 cells. Tracks measured in the movies of MCF10A cells were considered only if they were tracked for more than three frames (30 s) and tracks measured in movies of HL60 cells were only considered if they were tracked for more than three frames (6 s).

Statistical methods

Measurements of κ for MCF10A cells (Figure 4E) and actin-wave speeds for both cell types (Figure 5E) were compared on flat versus nanoridged surfaces using two-sample t tests with unequal variances. A two-tailed t distribution was used to calculate the reported p values. As the measurements κ of for HL60 cells violated the normality assumption for a t test, we used a nonparametric Wilcoxon rank sum test to compare these values. A full description of the statistical parameters involved in these tests is provided in Supplemental Dataset S1.

Data and software availability

Our optical-flow analysis code is available on GitHub at Imaging data from this study is available in a Mendeley Data repository on publication (Lee et al., 2020).


This article was published online ahead of print in MBoC in Press ( on February 5, 2020.

Abbreviations used:

extracellular matrix




multiphoton absorption polymerization




particle image velocimetry.


We thank the University of Maryland Imaging Incubator Core Facility for use of their systems in collecting images for this work. We appreciate Ema Smith’s work on the optical-flow code. This work was supported by AFOSR grant number FA9550-16-1-0052. R.M.L was supported by National Cancer Institute/National Institutes of Health (NIH) award number T32CA154274. L.C. was supported by COMBINE NRT award number 1632976 and NIH 1U01GM109887.


  • Azatov M, Sun X, Suberi A, Fourkas JT, Upadhyaya A (2017). Topography on a subcellular scale modulates cellular adhesions and actin stress fiber dynamics in tumor associated fibroblasts. Phys Biol 14, 065003. Crossref, MedlineGoogle Scholar
  • Barron JL, Fleet DJ, Beauchemin SS, Burkitt TA (1994). Performance of optical flow techniques. In: Proceedings 1992 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, New York: IEEE Computer Society Press, 236–242. Google Scholar
  • Begemann I, Saha T, Lamparter L, Rathmann I, Grill D, Golbach L, Rasch C, Keller U, Trappmann B, Matis M, et al. (2019). Mechanochemical self-organization determines search pattern in migratory cells. Nat Phys 1. Google Scholar
  • Blair D, Dufresne E (2018). Crocker-Grier Particle Tracking Algorithm. The Matlab Particle Tracking Code Repository. Available at: (accessed 9 August, 2018). Google Scholar
  • Bridges AA, Jentzsch MS, Oakes PW, Occhipinti P, Gladfelter AS (2016). Micron-scale plasma membrane curvature is recognized by the septin cytoskeleton. J Cell Biol 213, 23–32. Crossref, MedlineGoogle Scholar
  • Charras G, Sahai E (2014). Physical influences of the extracellular environment on cell migration. Nat Rev Mol Cell Biol 15, 813–824. Crossref, MedlineGoogle Scholar
  • Chen S, Hourwitz MJ, Campanello L, Fourkas JT, Losert W, Parent CA (2019). Actin cytoskeleton and focal adhesions regulate the biased migration of breast cancer cells on nanoscale asymmetric sawteeth. ACS Nano 13, 1454–1468. MedlineGoogle Scholar
  • Crocker JC, Grier DG (1996). Methods of digital video microscopy for colloidal studies. J Colloid Interface Sci 179, 298–310. CrossrefGoogle Scholar
  • Dalby MJ, Riehle MO, Sutherland DS, Agheli H, Curtis ASG (2005). Morphological and microarray analysis of human fibroblasts cultured on nanocolumns produced by colloidal lithography. Eur Cell Mater 9, 1–8 (discussion 8). Crossref, MedlineGoogle Scholar
  • Discher DE, Janmey P, Wang Y (2005). Tissue cells feel and respond to the stiffness of their substrate. Science 310, 1139–1143. Crossref, MedlineGoogle Scholar
  • Doyle AD, Wang FW, Matsumoto K, Yamada KM (2009). One-dimensional topography underlies three-dimensional fibrillar cell migration. J Cell Biol 184, 481–490. Crossref, MedlineGoogle Scholar
  • Driscoll MK, Sun X, Guven C, Fourkas JT, Losert W (2014). Cellular contact guidance through dynamic sensing of nanotopography. ACS Nano 8, 3546–3555. Crossref, MedlineGoogle Scholar
  • Graziano BR, Weiner OD (2014). Self-organization of protrusions and polarity during eukaryotic chemotaxis. Curr Opin Cell Biol 30, 60–67. Crossref, MedlineGoogle Scholar
  • Habermann B (2004). The BAR-domain family of proteins: a case of bending and binding? The membrane bending and GTPase-binding functions of proteins from the BAR-domain family. EMBO Rep 5, 250–255. Crossref, MedlineGoogle Scholar
  • Horn BKP, Schunck BG (1981). Determining optical flow. Artif Intell 17, 185–203. CrossrefGoogle Scholar
  • Inagaki N, Katsuno H (2017). Actin waves: origin of cell polarization and migration? Trends Cell Biol 27, 515–526. Crossref, MedlineGoogle Scholar
  • Ketchum CM, Sun X, Suberi A, Fourkas JT, Song W, Upadhyaya A (2018). Subcellular topography modulates actin dynamics and signaling in B-cells. Mol Biol Cell 29, 1732–1742. LinkGoogle Scholar
  • LaFratta CN, Fourkas JT, Baldacchini T, Farrer RA (2007). Multiphoton fabrication. Angew Chem Int Ed 46, 6238–6258. Crossref, MedlineGoogle Scholar
  • Lam Hui K, Kwak SI, Upadhyaya A (2014). Adhesion-dependent modulation of actin dynamics in Jurkat T cells. Cytoskeleton 71, 119–135. CrossrefGoogle Scholar
  • Lee R, Campanello L, Hourwitz M, Alvarez P, Omidvar A, Fourkas J, Losert W (2020). Data for: Quantifying topography-guided actin dynamics across scales using optical flow. Mendeley Data V1, doi: 10.17632/53ty3ryfdk.1. Google Scholar
  • Lee RM, Stuelten CH, Parent CA, Losert W (2016). Collective cell migration over long time scales reveals distinct phenotypes. Converg Sci Phys Oncol 2, 025001. Crossref, MedlineGoogle Scholar
  • Li L, Fourkas JT (2007). Multiphoton polymerization. Mater Today 10, 30–37. CrossrefGoogle Scholar
  • Lou H, Zhao W, Li X, Duan L, Powers A, Akamatsu M, Santoro F, McGuire AF, Cui Y, Drubin DG, et al. (2019). Membrane curvature underlies actin reorganization in response to nanoscale surface topography. Proc Natl Acad Sci 201910166. CrossrefGoogle Scholar
  • Lu P, Weaver VM, Werb Z (2012). The extracellular matrix: A dynamic niche in cancer progression. J Cell Biol 196, 395–406. Crossref, MedlineGoogle Scholar
  • Lucas BD, Kanade T (1981). An Iterative Image Registration Technique with an Application to Stereo Vision. In: Proceedings of the 7th International Joint Conference on Artificial Intelligence, Vol. 2, San Francisco, CA: Morgan Kaufmann Publishers, 674–679. Google Scholar
  • Maruthamuthu V, Aratyn-Schaus Y, Gardel ML (2010). Conserved F-actin dynamics and force transmission at cell adhesions. Curr Opin Cell Biol 22, 583–588. Crossref, MedlineGoogle Scholar
  • McMahon HT, Boucrot E (2015). Membrane curvature at a glance. J Cell Sci 128, 1065–1070. Crossref, MedlineGoogle Scholar
  • Melak M, Plessner M, Grosse R (2017). Actin visualization at a glance. J Cell Sci 130, 525–530. Crossref, MedlineGoogle Scholar
  • Meyer WH, Howard TH (1987). Actin polymerization and its relationship to locomotion and chemokinetic response in maturing human promyelocytic leukemia cells. Blood 70, 363–367. Crossref, MedlineGoogle Scholar
  • Miki H, Yamaguchi H, Suetsugu S, Takenawa T (2000). IRSp53 is an essential intermediate between Rac and WAVE in the regulation of membrane ruffling. Nature 408, 732. Crossref, MedlineGoogle Scholar
  • Müller J, Sixt M (2017). Cell Migration: Making the Waves. Curr Biol 27, R24–R25. Crossref, MedlineGoogle Scholar
  • Petrie RJ, Yamada KM (2015). Fibroblasts lead the way: a unified view of 3D cell motility. Trends Cell Biol 25, 666–674. Crossref, MedlineGoogle Scholar
  • Simoncelli E, Adelson EH, Heeger D (1991). Probability distributions of optical flow. In: Proceedings of the 1991 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, New York: IEEE Computer Society Press, 310–315. Google Scholar
  • Sun X, Driscoll MK, Guven C, Das S, Parent CA, Fourkas JT, Losert W (2015). Asymmetric nanotopography biases cytoskeletal dynamics and promotes unidirectional cell guidance. Proc Natl Acad Sci 112, 12557–12562. Crossref, MedlineGoogle Scholar
  • Sun X, Hourwitz MJ, Baker EM, Schmidt BUS, Losert W, Fourkas JT (2018). Replication of biocompatible, nanotopographic surfaces. Sci Rep 8, 564. Crossref, MedlineGoogle Scholar
  • Ventre M, Natale CF, Rianna C, Netti PA (2014). Topographic cell instructive patterns to control cell adhesion, polarization and migration. J R Soc Interface 11, 20140687. Crossref, MedlineGoogle Scholar
  • Vig DK, Hamby AE, Wolgemuth CW (2016). On the quantification of cellular velocity fields. Biophys J 110, 1469–1475. Crossref, MedlineGoogle Scholar
  • Vogel V, Sheetz M (2006). Local force and geometry sensing regulate cell functions. Nat Rev Mol Cell Biol 7, 265–275. Crossref, MedlineGoogle Scholar
  • Weiner OD, Marganski WA, Wu LF, Altschuler SJ, Kirschner MW (2007). An actin-based wave generator organizes cell motility. PLoS Biol 5, e221. Crossref, MedlineGoogle Scholar
  • Wolf K, Friedl P (2011). Extracellular matrix determinants of proteolytic and non-proteolytic cell migration. Trends Cell Biol 21, 736–744. Crossref, MedlineGoogle Scholar
  • Zimmerberg J, McLaughlin S (2004). Membrane curvature: how BAR domains bend bilayers. Curr Biol 14, R250–R252. Crossref, MedlineGoogle Scholar