Short-term bioelectric stimulation of collective cell migration in tissues reprograms long-term supracellular dynamics

Short-term bioelectric stimulation of collective cell migration in tissues reprograms long-term supracellular dynamics

Abstract

The ability to program collective cell migration can allow us to control critical multicellular processes in development, regenerative medicine, and invasive disease. However, while various technologies exist to make individual cells migrate, translating these tools to control myriad, collectively interacting cells within a single tissue poses many challenges. For instance, do cells within the same tissue interpret a global migration ‘command’ differently based on where they are in the tissue? Similarly, since no stimulus is permanent, what are the long-term effects of transient commands on collective cell dynamics? object_of_studyWe investigate these questions by bioelectrically programming large epithelial tissues to globally migrate ‘rightward’ via electrotaxis. sensescontent_tagsTissues clearly developed distinct rear, middle, side, and front responses to a single global migration stimulus. claim_typephenomenonFurthermore, at no point poststimulation did tissues return to their prestimulation behavior, instead equilibrating to a 3rd, new migratory state. These unique dynamics suggested that programmed migration resets tissue mechanical state, which was confirmed by transient chemical disruption of cell–cell junctions, analysis of strain wave propagation patterns, and quantification of cellular crowd dynamics. contributionOverall, this work demonstrates how externally driving the collective migration of a tissue can reprogram baseline cell–cell interactions and collective dynamics, even well beyond the end of the global migratory cue, and emphasizes the importance of considering the supracellular context of tissues and other collectives when attempting to program crowd behaviors.

Coordinated cell migration is critical in processes spanning development, healing, and disease; accordingly, improving our ability to program and engineer such collective migration is an important goal. A key challenge is to understand how a group responds to given migration ‘commands’ as a collective unit, particularly as such cues are often transient. We investigated this by stimulating epithelial tissues to migrate ‘rightward’ using a bioelectric stimulus, electrotaxis, and mapped how different parts of the tissue responded to the same command, both during and after stimulation. We find that cells within stimulated tissues respond based on where they are in the tissue, and that such global commands induce reprogramming of innate collective cell behaviors even long after the stimulus stops.

Introduction

Collective cell migration is essential for many multicellular organisms and underpins diverse processes, spanning organ development to cancer invasion ( 1–3 ). Therefore, externally directing collective cell migration should allow us to not only better understand these processes, but also to formally program and engineer them for applications in regenerative medicine and tissue engineering, where large scale cellular motion is necessary, such as in wound healing ( 4–9 ). concepts↗concepts↗content_tags↗Programming cell migration is increasingly feasible with myriad emerging approaches spanning chemotactic devices (cells migrate along chemical gradients) ( 8 , 10–14 ), light-induced, optogenetic modulation of signaling pathways ( 15–17 ), and bioelectric programming of cell migration through electrotaxis (cells migrate along electrical gradients) ( 9 , 18–24 ). However, a key challenge for any stimulation modality is understanding how global commands, such as ‘All Cells Migrate Rightward’, play out and drive cellular dynamics at the overall system or tissue level where 10,000+ cellular agents must coordinate at the macroscale.

It is well established that cell–cell interactions can cause a different stimulus response to that observed with single cells. For instance, studies have demonstrated both that cellular groups often better follow chemical ( 12 ) and electrical migration cues ( 25 ) than do single cells, while overly strong cell–cell interactions and cell–cell coupling can actually compete with external migratory commands leading to adverse outcomes ( 7 , 26 ). research_questionconcepts↗Motivated by these considerations, we investigated 2 key factors for externally directing collective cell migration: (1) how does where a cell is located within a tissue modulate its response to a global command; (2) does following this global command change the longer-term behavior of the tissue after the command is removed?

It is increasingly clear across both in vitro and in vivo studies of collective cell migration that the location of a cell within a group strongly affects how that cell behaves ( 27–33 ). In the tissue context, a growing body of work has revealed how gradients of mechanical tension underlie the growth and motion of highly collective tissues, such as the epithelia lining our organs ( 34–37 ), and naturally give rise to behavioral differences between boundaries, edge zones, and bulk of a tissue ( 27 , 29 ). Such behavioral zones can emerge even in groups of loosely coupled cells migrating along chemical gradients ( 38 ). concepts↗This compartmentalization of behaviors within a group is referred to as ‘supracellularity’ ( 38 ), and emphasizes how collectives are often best characterized as a whole unit. Applying this supracellular framework here provides an analytical basis to understand how, when, and where cells within a single tissue respond to the same global migratory command.

Additionally, since no migratory command or stimulus is permanent, it is critical to understand the collective response of a cellular group or tissue not only during a global migratory stimulation, but also after removal of the stimulus. Considering the whole process, from initial entrainment of cells to eventual relaxation, can reveal the longer-term consequences of externally driving collective cell behaviors, and can detect a collective ‘memory’ of the stimulation event. To date, such cellular memory of migratory stimuli has been studied primarily with single cell relaxation after exposure to chemotactic or electrotactic cues ( 10 , 20 ). After a chemotactic gradient is removed, internal signaling such as Ras activity may relax over a 30-second scale ( 39 ), while front-rear cell polarity and overall directionality decay within ∼10 minutes poststimulation ( 10 , 11 ). In electrotaxing single cells, directionality has been reported to persist over a ∼15–60-minute period poststimulation ( 20 , 40 ). However, the collective nature of a tissue relative to single cell studies means tissues will necessarily exhibit different relaxation and reprogramming dynamics in response to a global migration command; investigating the spatiotemporal response of the tissue to external commands will go a long way to improving both our understanding of collective cell behaviors, and our ability to more effectively program and augment large-scale tissue behaviors for tissue engineering.

To investigate these questions, we needed both a model system and an appropriate migration stimulus. concepts↗We began with a gold-standard collective cell migration model—the MDCK renal epithelium, which is one of the most well-characterized collective cell systems to date ( 41–43 ). MDCK epithelia are readily grown to centimeter scale, and exhibit robust outward migration and canonical ‘scratch wound’ healing dynamics ( 27 , 29 , 44 ), enabled by strong cell–cell coupling mediated by E-cadherin and other junctional complexes ( 18 , 35 , 45 ). We next needed a compatible stimulus that could precisely and globally induce directional migration across the whole epithelium. Chemotaxis and electrotaxis are 2 of the best characterized directional migration candidates. However, while chemotaxis is potent and well-understood in certain single cells and small cluster models, many cell types and model systems lack chemotactic pathways or receptors (including MDCK cells), so chemotaxis needs to be engineered and fine-tuned in cells that do not natively chemotax ( 8 ). Electrotaxis, by contrast, arises naturally in most cell types and is conserved across many species ( 20 , 24 , 46–48 ). Electrotaxis shares certain key signaling pathways with chemotaxis (notably PI3K and TORC2) ( 48 ), and induces directional migration based on naturally occurring ionic current gradients ( \({\mathcal O}\) (1 V/cm)) ( 49 , 50 ) formed during tissue perturbations (e.g. morphogenesis and healing). These fields apply weak electrophoretic and electro-osmotic forces to membrane receptors, inducing receptor aggregation and activation, leading to cell polarity and migration ( 20 , 51 ). That electrotaxis affords programmatic control of collective cell migration with uniform and near-instantaneous stimulus of the whole tissue ( 19 , 26 , 52 , 53 ), without requiring additional reagents or cellular modifications, made it an optimal ‘command’ stimulus for this study.

We then used automated phase contrast microscopy to capture supracellular patterns and relaxation dynamics of epithelial tissues exposed to a transient electric field inducing ‘rightward’ electrotaxis. We first quantify how the various edges (leading, trailing, top, and bottom) of a ‘rightward’ migrating tissue exhibit distinct behaviors, with particular emphasis on apparent elastic recoil in these zones once stimulation is removed. These responses were markedly different from how the bulk center of tissues behave, where we observed long-term directional ‘memory’ long after stimulation was stopped. We link this memory to a resetting of tissue mechanical state that must occur in order to enable collective migration of mechanically coupled cells, which we validate by characterizing strain wave dynamics before, during, and after stimulation. Finally, we demonstrate how driven collective migration ultimately reprograms the underlying cell–cell correlations in a lasting manner, meaning that a poststimulation tissue differs markedly from its unstimulated, control sample counterpart.

Results

Programmed collective migration alters supracellular migration patterns in tissues

To generate the core data needed to explore supracellularity and programmed migration, we captured the complete step response of large epithelia to global bioelectric migration commands (Fig. 1A – C , Movie 1 compares control and stimulated tissues). method↗This required continuous time-lapse imaging of arrays of 5 × 5 mm MDCK epithelial monolayers over 3 contiguous periods: control (1 h), stimulation ON (3 h), and stimulation OFF (6 h). Movie 2 shows closer views of local cellular responses during this process in control and stimulated tissues. To induce electrotaxis, we built custom electro-bioreactors around these tissue arrays similar to our SCHEEPDOG platform ( 7 , 19 ), which allowed us to provide continuous media perfusion and a computer-controlled, uniform, ‘rightward’ electrical cue (3 V/cm field) to all tissues in the bioreactor. A schematic of our bioreactor and representative epithelium highlighting our analysis zones are presented in Fig. 1A and B .

concepts↗concepts↗concept_edgesTo begin, we measured collective migration patterns using particle image velocimetry (PIV) on phase contrast images (PIV; see Methods) to map the spatial velocity profile across the tissue and throughout the experiment (Fig. 1C ). We report velocities parallel ( V x , Fig. 1D – F ) and perpendicular ( V y , Fig. 1G – I ) to the direction of stimulation (full heatmaps presented in Movie 3 and Figure S1, Supplementary Material ). In the direction parallel to the field, we saw clear progression of the velocity field, from relatively disordered during the control period (Fig. 1D ), to higher magnitude and well-aligned with the field during stimulation (Fig. 1E ), and finally to a distinct relaxation phase poststimulation (Fig. 1F ). In the direction perpendicular to the field (Fig. 1G – I ), we observed no obvious change during stimulation, but a similar poststimulation state to Fig. 1F , except parallel to the y -axis.

Collective electrotaxis step response induces a spatiotemporal effect in epithelial tissues persisting long after stimulation. ( A ) Diagram of electro-bioreactor providing electrical stimulation in the x -axis and perfusion of fresh media in the y -axis through the chamber containing three 5 × 5 mm tissues. ( B ) Phase contrast image of a single MDCK epithelial tissue; annotations denote specific regions analyzed throughout the paper. ( C ) Trace of electrical stimulation step response stimulus: 1 h unstimulated, 3 h ‘rightward’ stimulation, and 6 h poststimulation. ( D )–( F ) Single tissue heatmaps of velocity in the x -direction, aligned with the direction of electric field (red). Note that (E) indicates large increase in the rightward collective motion during stimulation, while (F) indicates a new relaxation behavior hours after stimulation. See Figure S1 (Supplementary Material) . for corresponding control tissue heatmaps. ( G )–( I ) Analogous heatmaps of velocity in the y -axis, perpendicular to the direction of electric field. Note again the new pattern long after stimulation ended in (I), despite minimal apparent changes during stimulation in (H).

While directed collective migration during stimulation is expected with tissue electrotaxis ( 4 , 18 , 19 , 52 ), the distinct migratory patterns long after stimulation had ceased, and the fact that these patterns did not resemble any previous state were surprising and previously unreported. Furthermore, the migration patterns visible in poststimulation heatmaps, where the bulk zone retained a ‘memory’ of rightward migration (red center, Fig. 1E and F ) while the perpendicular and trailing edges appeared to ‘recoil’ and migrate opposite to the prior command (blue edges, Fig. 1F ), suggest a long-term reprogramming of the supracellular dynamics. These compartmentalized behavioral zones indicated where to focus our further analyses.

Tissue boundaries compete with external migration commands and rapidly recover poststimulation

We separately analyzed the different edge and bulk zones to better understand how a tissue parses a global migration command, beginning with the leading and trailing edge dynamics, as edge outgrowth is very well-characterized in unperturbed epithelia ( 27 , 29 , 42 ). content_tagsHere, we found that while both the leading and trailing edges obeyed the migration command and exhibited greater ‘rightward’ migration speed during stimulation, they very rapidly equilibrated back to control growth speeds in < 30 minutes poststimulation (Fig. 2A ; see also Fig. 1D – F ). Interestingly, the trailing edge exhibited a transient overshoot of outgrowth speed that was quickly damped within 30 minutes poststimulation, suggestive of the viscoelastic nature of epithelial tissues ( 54 , 55 ). To resolve how individual cellular agents behaved in the edge zone, we performed single-cell tracking of all cells within 150 μm of the trailing edge and overlaid cell trajectories in a ‘hairball’ plot (Fig. 2B ; see Methods). This clearly demonstrates both directed migration following the stimulus (pink line segment of Fig. 2B ), and how cells performed rapid ‘U-turns’ immediately poststimulation and quickly equilibrated to control edge dynamics.

Edges parallel and perpendicular to the field direction displayed distinct behaviors. Cartoons indicate analysis zones. ( A ) Edge outgrowth speed (positive is rightward motion). Control edges (gray) exhibit steady outwards motion, while stimulated edges exhibit perturbed outgrowth during stimulation with rapid poststimulation recovery. ( B ) Cell trajectories from the trailing zone. Lines color-coded by time (see color bar). Mean trajectory of all tracked trailing zone cells indicated with bold gray line; magenta overlay indicates period of active stimulation. All trajectories were plotted as positive final y -displacements. See Movie 3. ( C ) Average speed in the direction of the electric field in top and bottom edge zones. Note strong increase in rightward motion in these zones during stimulation and rapid over-correction poststimulation (antiparallel motion). ***: P = 0.0002 between stimulated and control tissues over 6–10 h. ( D ) Cell trajectories from top edge zone. Note again recoiling motion poststimulation. For panels (A) and (C), shading indicates SD across N control = 6 and N stimulated = 9; panels (B) and (D) represent single tissues, where zones exclude 1.5 mm on each side for edge effects.

More surprising was the behavior of the top and bottom edges of the tissue. As indicated in Fig. 1 , these edge zones largely followed the electrotactic command during stimulation (Fig. 1E and H ; Figure S1, Supplementary Material) , but then exhibited retrograde motion poststimulation as the epithelium re-established conventional expansion into the free surrounding space, qualitatively reminiscent of elastic recoil. Again, analyzing a 150-μm zone at the top and bottom edges, we plotted the ensemble migration velocity (from PIV) in the direction of the stimulus (Fig. 2C ; see also Fig. 1 and Figure S1, Supplementary Material ). While cells in both top and bottom border zones migrated in the programmed ‘rightward’ direction and increased speed during stimulation, their poststimulation dynamics exhibited strong, persistent ‘leftward’ migration over several hours (Fig. 2C ; Figure S2, Supplementary Material ; P = 0.0002 vs. control over 6–10 h; see Methods). A hairball plot of cells in the ‘top’ zone (Fig. 2D ) again demonstrates initial entrainment to the electrical command (pink line segment of Fig. 2D ) followed by a distinct ‘recoil’ poststimulation as the tissue re-equilibrates to the conventional outward migration. The dynamics of this process are highlighted in Movie 4, showing a tracked group of cells near the top boundary (see Methods). mechanisms↗Overall, these distinct edge dynamics demonstrate that external commands compete with the strong outward migration program of epithelial tissue edges, resulting in a local tug-of-war, eventually won by the outward migration program poststimulation. Further, this apparent recoil in the nonleading edges is suggestive of elastic recovery after removal of a stress caused by electrical stimulation, consistent with prior work indicating that electrotaxis notably elevates intercellular stresses in the bulk zone of a tissue ( 56 ).

Cells in the tissue center exhibit a memory of the prior migratory command

sensesIn contrast to the distinct edge recoil behaviors, the bulk of the tissue relaxed far more gradually and appeared to be the primary source of migratory ‘memory’. To best emphasize this, we analyzed the bulk zone using standard ensemble analyses for directed collective cell migration—the directionality order parameter and mean speed ( 18 , 25 , 26 , 57 ), both plotted over the entire experimental time course. Here, we define the directionality order parameter as (1) \[\begin{eqnarray*} \varphi = \frac{1}{n} \sum cos(\theta _i), \end{eqnarray*}\] where θ is the angle between a PIV velocity vector and the direction of stimulation, the ‘command’ direction unit vector. content_tags↗φ ranges from [−1,1], with −1 and 1 indicating cells moving exactly antiparallel and parallel to the command direction, respectively. φ = 0 indicates no net directionality at the ensemble level, as is the case for the control data shown (Fig. 3A , black dashed line). In stimulated tissues (Fig. 3A , red line), the directionality dynamics clearly resolve the fast transition from control ( φ ∼ 0) to stimulated ( φ ∼ 1), followed by a much slower relaxation of the tissue migration direction poststimulation. In contrast, while speed in the direction of stimulation, | V x |, increased nearly 5× during electrotaxis (Fig. 3B ), this effect rapidly decayed to a level that was actually lower than the speed of an equivalent control tissue (Fig. 3B , P = 0.003 for 6–10 h; see Methods).

Tissue centers show migratory memory poststimulation. All data here are captured from the tissue bulk region (see cartoon). ( A ) Directionality order parameter φ comparison between unstimulated (dashed gray) and stimulated (red) tissues. Note sustained directionality even poststimulation. **: P = 0.003 between mean values of stimulated and control tissues over 6–10 h. ( B ) Average horizontal speed (parallel to field) during stimulation. Note both strong rise during stimulation and sharp fall poststimulation to below control levels. **: P = 0.003 over 6–10 h. ( C ) and ( D ) Exponential decay modeling of poststimulation relaxation dynamics τ for low, medium and high density tissues (see Methods). Raw data in (C) and (D) indicated by dots; solid curves represent best fits. ( E ) Quantification of characteristic timescales with respect to cell densities for directionality (left) and horizontal speed (right). Note slow directionality relaxation and inverse relation between density and directionality relaxation, which is not true for speed relaxation. For panels (A) and (B), N control = 6 and N stimulated = 9; shading indicates SD; (C)–(E) include N = 5, 9, and 4 for low, medium, and high density tissues, respectively. In all panels, tissues’ center 2 × 2 mm were analyzed.

To quantify this difference in relaxation times between ensemble directionality and speed, we defined these timescales τ by fitting the data in the poststimulation phase to an exponential decay model (Fig. 3C and D ; see Methods). concepts↗content_tagsHere, the characteristic relaxation time τ for the directionality order parameter (∼160 min) was greater than 5-fold longer than for speed (∼30 min), as shown in Fig. 3C – E (red). To put this in context, these epithelial cells lost ∼80% of their previous speed within 20 minutes (losing ∼1 μm/min), yet maintained an imperative and ‘memory’ to continue migrating ‘rightward’ for many hours (until φ relaxed closer to 0). To add further biophysical context, and motivated by the role cell density plays in modulating collective cell migration and tissue mechanics ( 29 , 58 , 59 ), we repeated these relaxation studies in the tissue bulk at 2 additional densities, the lower bound being the minimum density that maintained confluence, ∼2,000 cells/mm 2 (Fig. 3C – E , orange), and the upper bound, ∼4,500 cells/mm 2 (Fig. 3C – E , brown), sufficient to geometrically limit migration (see Methods) ( 60 ). Relaxation dynamics versus density are shown in Fig. 3C – E . content_tagsWhile speed relaxed independently from density, migration directionality relaxation times were inversely correlated with density (Fig. 3E ), ranging from ∼3 h (Fig. 3E , low density) to ∼1 h (Fig. 3E , high density). We hypothesized that this trend may be due to alterations in cell shape and overcrowding at higher densities ( 3 , 59–63 ).

To specifically investigate if global electrotaxis might reprogram cell shape and, in turn, cell shape modifications might alter the system response, we analyzed both the cell shape index and cellular orientation shifts. concepts↗Cell shape index is a standard geometric benchmark for assessing cellular jamming and tissue mechanical state ( 62 , 64 , 65 ), defined as \(q = \left\langle Perimeter_i/\sqrt{Area_i}\right\rangle\) , where the brackets denote an average taken over this ratio for all cells. content_tags↗content_tagsGenerally, an average cell shape index above the critical value of 3.81 in epithelia indicates a more fluid-like behavior, while a shape index below 3.81 indicates a more solid-like tissue. Interestingly, when we assessed the average shape index in the bulk of electrotaxing tissues using a modified Voronoi approach ( Figure S3, Supplementary Material ; see Methods) ( 66 , 67 ), concepts↗content_tagswe found no indication of a systemic jamming transition induced by stimulation ( Figure S4, Supplementary Material ), meaning that neither the eventual reductions in speed nor relaxation and memory dynamics were explicitly the result of a jamming transition, even in our densest cases. It should also briefly be noted that cell death was not a contributing factor here, since total cell count showed no sign of cell death in stimulated tissues compared to controls ( Figure S5, Supplementary Material ). That said, there was a notable decrease in shape index following stimulation that exceeded the expected decrease in the control case ( Figure S4, Supplementary Material ; poststimulation). limitationsWhile this is perhaps suggestive of a mechanical recovery occurring poststimulation resulting in a slight increase in solid-like properties in the tissue, future work is necessary to further elucidate this.

To better capture how the shape index behaved spatiotemporally throughout the tissue, rather than simply in the bulk, we created kymographs of cell shape index over the entire tissue width for 3 different densities ( Figure S6, Supplementary Material ; see Methods). While expected variations were present at the edges of tissues as cells spread into free space, we observed no obvious spatial patterning resulting from electrotaxis, beyond observing the same apparent decrease in shape index poststimulation.

Cell shape index does not account for potential anisotropy in cellular shape and potential alignment, and a hallmark of electrotaxis in many cell types is that the long axis of cells tends to align perpendicular to the field axis ( 56 ). We investigated that here by using a version of our MDCK-II cells that stably expressed E-cadherin:DsRed ( 68 ) to determine both if cell shape alignment occurred in our tissues during stimulation, and if there was a poststimulation memory of shape anisotropy. Fluorescent E-cadherin is an excellent boundary marker in epithelial cells, and allowed us to observe electrotaxis behavior at higher resolution (20× vs. 5× previously) to capture subtleties in membrane shape and electrotactic dynamics (Movie 5 highlights this). First, we analyzed cell orientation using our labeled cells and orientation mapping (see Methods) and observed the expected long-axis alignment perpendicular to the field ( Figure S7, Supplementary Material ; middle column) ( 56 ). While this realignment occurred in both confluent (Movie 5, bulk and edge zones) and single cells (Movie 5, bottom), only larger clusters exhibited directional migration, supporting prior studies indicating that epithelial electrotaxis may be an emergent phenomenon at the collective level ( 26 ). More surprising was the preliminary analysis of cell orientation poststimulation. Here, the poststimulation data in the tissue bulk ( Figure S7, Supplementary Material ; top right) still indicate a biased long axis, suggesting that electrical stimulation induced a long-term realignment of cells in the bulk. This trend was not obvious in cells either at the edges of a monolayer or in small clusters, indicating both that neighbors are needed for stronger orientation memory and again that boundary dynamics differ from bulk dynamics. Specifically, the slow relaxation of cell orientation and directionality in the bulk (Fig. 3 ) is quite distinct from the rapid poststimulation response in the tissue edges (Fig. 2 ), perhaps suggesting a decoupling of the bulk from the edges. To investigate this further, we next characterized interplay between edge and bulk dynamics.

Halting bioelectric stimulation resets the mechanical state of the tissue

To better capture the supracellular dynamics’ coupling behaviors from the edges to the bulk, we generated average kymographs of the directionality order parameter, φ , across the entire width of the tissue (averaged vertically; see Methods) over the entire experimental duration (Fig. 4A – C , red again indicates rightward motion). We first established a control case as a reference, Fig. 4A , demonstrating canonical, steady outward tissue expansion in both directions. By contrast, while electrical stimulation induced an overwhelmingly ordered response during the stimulation period (Fig. 4B , red center in 1–4 h; see also Fig. 1E ), the poststimulation dynamics revealed pronounced, inward traveling waves of migration mobilization that nucleated at the outer edges of the tissue and propagated inward at an apparently constant rate (‘triangular regions’ marked with dotted lines in Fig. 4B ). Such waves also clearly emerged in poststimulated tissues from the top and bottom edges in V y kymographs ( Figure S8, Supplementary Material ). Generally, inward traveling epithelial waves have previously been reported as a natural consequence of releasing a confluent epithelium from confinement and allowing it to grow outward (as in Fig. 4A ) ( 37 , 57 , 69 , 70 ). However, the poststimulation effect observed here was far stronger than what we observed in our control tissues (compare Fig. 4A and B ).

Kymographs reveal apparent mechanical reset postelectrical stimulation. ( A )–( C ) Averaged directionality order parameter φ kymographs (see Methods) demonstrate the evolution of x -axis motion patterns across a whole tissue. (A) Control tissues, (B) electrically stimulated tissues, and (C) tissues with EGTA-disrupted cell–cell junctions (see Methods), where the end of treatment was timed to coincide with the end of stimulation in (B) for optimal comparison. Dashes diagonal lines in (B) and (C) demarcate the ‘triangles’—inward traveling waves of outward migration induced by mechanical reset. ( D ) Plot of the positions of these traveling waves in (B) and (C); lines are linear best fits, see panel legend and Methods. ( E ) Traveling wave speeds from slopes of linear fits in (D): EGTA junction disruption (green); and electrically stimulated tissues (gray). Error bars indicate SD. ( F ) Profile at t = 7 h of directionality order parameter φ across all 3 kymographs, as marked by the dashed black line across panels (A)–(C). Shading indicates SD across tissues. For all panels, N control = 6, N stimulated = 9, and N EGTA = 3.

We hypothesized that these large-scale inward traveling waves indicated that removal of the global electrical cue significantly altered the mechanical state of the tissue. It has previously been demonstrated that steady-state electrotaxis can alter tissue mechanics with increased traction force magnitude in a tissue center and reoriented traction force alignment perpendicular to the electric field vector ( 56 ), hence, sudden removal of the electrotactic cue should also alter mechanics. However, a better analog to the large, inward waves we observed comes from a prior study without electrical stimulation that demonstrated that transient disruption of cell-cell adhesion within an epithelium appeared to reset the internal mechanical tension and cell traction state of the tissue, while also inducing similar large-scale inward waves ( 37 ). For comparison, we replicated this perturbation using a 30-minute pulse of EGTA (Fig. 4C ; see Methods) to briefly chelate calcium and transiently disassemble E-cadherin binding between cells in an otherwise control tissue, and observed remarkably similar inward traveling waves of cell mobilization upon restoration of calcium. mechanisms↗content_tagscontent_tagsWe confirmed that the inward traveling waves induced by both electrical stimulation and calcium chelation featured linear propagation at rates of ∼100 μm/h (Fig. 4D and E ), in agreement with those previously described ( 37 ) and consistent with a large-scale mechanical ‘reset’ caused by removal of the global bioelectric command.

While halting electrotaxis and transiently disrupting cell–cell adhesion produced similar effects on tissue boundaries, a key difference in the response lies in the bulk tissue behavior, where poststimulation tissues still exhibited directional memory in the bulk (Fig. 3 , center of 4B in 4–10 h range), while chemical disruption of cell–cell junctions largely shut down cell motility in the bulk (center of Fig. 4C in 4–10 h range). We quantified this by analyzing line sections taken across each kymograph at t = 7 h (Fig. 4F ), giving us a snapshot of the profile across the tissues 3 h after the respective perturbations. Control tissues exhibited a smooth and largely symmetric, graded directionality profile from the center of the tissue outward (Fig. 4F , black dashed curve). Tissues with transient chemical disruption of cell–cell junctions exhibited highly directed edge zones, but no net directionality in the bulk (Fig. 4F , green curve). However, poststimulation tissues exhibited not only similar directed edge zones to the junctional disruption case, but a central zone with net positive directionality, implying persistent rightward motion long after stimulation had ceased (Fig. 4F , red curve, matching ‘memory’ shown in Fig. 3 ). Hence, ending electrical stimulation appears to reset the tissue’s mechanical state at the boundaries similarly to chemically disrupted cell–cell adhesion. However, despite similarity in migratory behavior at tissue edges, it is unlikely that electrotaxis disrupts cell–cell junctions as E-cadherin has been shown to be necessary for MDCK and certain other epithelial electrotaxis ( 25 , 26 ). content_tagsMoreover, the migratory memory and preservation of front–rear polarity in the tissue bulk was only observed postelectrotaxis, and not postjunctional disruption.

Programmed migration with electrotaxis disrupts collective strain wave propagation

As commanding a tissue to migrate ‘rightward’, and abruptly negating that command, both appear to alter the tissue mechanical state in a location dependent fashion akin to junctional disruption, we next sought a biophysical mechanism coupling collective migration behaviors across a tissue. Critically, the waves from Fig. 4 represent coordinated inward-traveling mechanical waves of outward-directed migration, which ultimately allow the motile tissue edges to mechanically influence migration within the tissue bulk. For instance, if a leading edge cell in a cohesive tissue migrates outward, then its immediate rearward-adhered neighbor can follow, allowing that cell’s own rearward neighbor to follow suit, and so on, thus forming a wave of force-coupling and cell mobilization. concepts↗Such waves have been shown to create zones of stretching and compressing, that are transmitted across the tissue via ‘strain waves’, and tend to originate at tissue boundaries ( 37 , 59 ). To capture these waves and analyze the overall rate of cellular deformations, we measured the rate of strain—stretching (positive) or compressing (negative)—from the PIV vector fields, as \(\frac{dV_x}{dx}\) for strain rate in the x -direction ( \(\dot{\epsilon }_{xx}\) ) and \(\frac{dV_y}{dy}\) for strain rate in the y -direction ( \(\dot{\epsilon }_{yy}\) ).

Since strain waves are difficult to visualize in such large tissues, we represented our control tissue bulk data as strain rate kymographs, a horizontal strip for \(\dot{\epsilon }_{xx}\) (Fig. 5A ) and a vertical strip for \(\dot{\epsilon }_{yy}\) (Fig. 5B ), from the strain rate vector fields, similar in form and presentation to the strain waves depicted in prior studies ( 37 ). Waves of \(\dot{\epsilon }_{xx}\) traveled primarily in the x -direction, appearing as diagonal regions of stretch (tensile strain rate; purple in Fig. 5A ) or compression (compressive strain rate; green in Fig. 5A ) in the x–t kymograph, with the slopes confirming that these waves propagated at similar rates as the triangular waves previously discussed in Fig. 4 . Similarly, waves traveled primarily in the y -direction, seen in the y–t kymograph (Fig. 5B , also depicted with sample slope speed illustrations). Movie 6, with control and stimulated tissues side-by-side, visually portrays general changes in strain rate dynamics, especially the boundary and bulk effects poststimulation.

Strain waves reorient based on electrical stimulation. ( A ) Representative x–t kymograph (1 mm over 10 h) of bulk strain rate \(\dot{\epsilon }_{xx}\) in unstimulated tissues showing horizontal strain wave propagation. Positive strain rate (purple) indicates stretch; negative strain (green) indicates compression. Wave propagation speeds (slope of dashed orange lines) varied from 75 to 120 μm/h. ( B ) Representative y–t kymograph of bulk strain rate \(\dot{\epsilon }_{yy}\) ; note vertical orientation relative to 5A to emphasize vertical strain wave propagation. ( C ) Polar histograms of strain rate wave direction, θ xx , in stimulated tissues plotted for 3 experimental regimes: unstimulated (gray, 0–1 h); stimulated (red, 2–4 h); and poststimulation (yellow, 8–10 h). Note sustained leftward shift of θ xx poststimulation. ( D ) Polar histograms of strain rate wave direction, θ yy , analogous to 5C . Note transient rightward orientation during stimulation (red), followed by recovery to baseline poststimulation (yellow) strongly matching the unstimulated case (gray). ( E ) Quantification of average strain wave directionality (mean cos ( θ ii ); −1 is leftward, +1 is rightward) for each condition. Control tissue cases of \(\dot{\epsilon }_{xx}\) (black) and (gray) show no net directionality. Orange squares show mean cos ( θ xx ) time-course in stimulated tissues. Note perturbation during stimulation and persistent lack of recovery, biased leftward. Purple triangles show mean cos ( θ yy ) in stimulated tissues. Note rightward perturbation during stimulation and clear recovery poststimulation. Shading indicates SD across tissues. Panels (C)–(E) data is averaged over N = 9 stimulated tissues.

Kymographs of strain rate only visually capture wave propagation that persists for several hours. To quantify the propagation direction of these waves in the shorter timescales of our experiments, we performed additional PIV analysis on image sequences of strain waves themselves (see Methods). This produced vector fields of the primary propagation direction of strain waves at each location within the tissue for every timepoint, which we denote as θ xx (for \(\dot{\epsilon }_{xx}\) waves) and θ yy (for waves). We then plotted polar histograms of θ xx and θ yy within the tissue bulk to visualize the distribution of strain rate propagation direction (Fig. 5C and D ). In the unstimulated control hour (0 – 1 h; gray-shaded left histograms of Fig. 5C and D ), as in control tissues (above, Fig. 5A and B ), \(\dot{\epsilon }_{xx}\) waves travel primarily horizontally with equal leftward- and rightward-moving waves (Fig. 5C , left), and waves travel primarily vertically with balanced upward- and downward-moving waves (Fig. 5D , left). This matches expectations, whereby information in a symmetric, unperturbed tissue should propagate equally from all free edges throughout the tissue.

Surprisingly however, these dynamics dramatically changed during and after electrical stimulation, where the \(\dot{\epsilon }_{xx}\) strain waves primarily propagated leftward, with this effect being even more pronounced poststimulation than during electrotaxis itself (Fig. 5C , center and right). This means that globally ‘rightward’ electrotaxis induces a mechanical tissue state, such that horizontally traveling waves within the tissues are dominated by those that would normally stem from the rightward (leading) edge. By contrast, \(\dot{\epsilon }_{yy}\) waves were strongly biased in the direction of stimulation during electrotaxis (90° of their baseline orientation), essentially shifting with the cell migration itself, but this effect relaxed poststimulation to match control orientation (Fig. 5D , center and right). These data suggest that, while active electrotaxis appears to reprogram much of the endogenous mechanical strain state within a tissue, the lasting impact lies only along the axis of induced migration.

To further illuminate this, we analyzed the dynamics of the strain wave disruption process using the average value of cos ( θ xx ) and of cos ( θ yy ) within the tissue bulk as a metric for the directionality of \(\dot{\epsilon }_{xx}\) and waves, respectively (Fig. 5E ). A directionality of 1 would represent a tissue for which waves are only traveling parallel to the field, while a directionality of −1 would represent waves only traveling antiparallel. \(\dot{\epsilon }_{yy}\) directionality peaked to the right about 1 h after stimulation was initiated and decayed quickly after stimulation was turned off (Fig. 5E , purple triangles), a time-course which is strikingly reminiscent to that of speed (Fig. 3B ). \(\dot{\epsilon }_{xx}\) waves were more nuanced, with large negative jumps in directionality in the first timepoint of stimulation and the first timepoint poststimulation, each followed by dynamic changes which are perhaps tied to viscoelastic effects. Poststimulation, however, the \(\dot{\epsilon }_{xx}\) directionality eventually reaches a distinct new steady state—a long-term bias toward the trailing edge with no apparent decay (Fig. 5E , orange squares). The initial immediate changes to wave propagation bear resemblance to waves of signaling that span a tissue in the matter of minutes as a response to tissue damage, which also promptly move leftward from a wound at the right edge ( 71–73 ), in contrast to the much slower mechanical waves that produced the boundary ‘triangles’ in the kymographs in Fig. 4B and C and of migration upon release of a barrier ( 37 ). In these data, we see that both initiation and removal of stimulation strongly reprogrammed the mechanical state of the tissue. This long-lasting change of mechanical communication should produce changes in the collective structure of cell migration, which we analyze next.

Driving collective cell migration overrides and weakens cell – cell crowd interaction dynamics even poststimulation

A key hallmark of any collective motion process is that information couples well beyond simple nearest-neighbor interactions. In bird flocks, neighbor correlation dynamics may propagate up to 7 neighbors ( 74 ), while the MDCK epithelium is known to exhibit correlated domains of 5 – 10+ cells ( 28 , 75 ), mediated by cell – cell adhesion ( 76–78 ). However, programming a particular pattern of migration into a group that already possesses internal collective coupling necessitates a conflict ( 7 ), mediated in electrotaxing epithelia by E-cadherin cell – cell adhesion ( 4 , 7 , 25 ). Given this framework, we hypothesized that electrically programming collective cell migration might reprogram the endogenous correlation length to force cells to entrain to the migration command, and that the relaxation behaviors we observed were driven by a re-establishment of the native correlation dynamics.

To test this hypothesis, we calculated the velocity – velocity spatial correlation length both parallel ( V x ) and perpendicular ( V y ) to the electric field stimulation axis at each time point using PIV data. concepts↗content_tags↗The correlation length reflects the approximate size of correlated domains within a larger population. When calculating velocity – velocity correlations in a highly directed system (e.g. flocks of birds, road traffic, or ensembles of electrotaxing cells), it is critical to first subtract the global mean velocity at any given time, and therefore, perform correlation analysis on the velocity residuals, or fluctuations about the mean, which much better capture the internal cell – cell interactions by removing the global bias ( 42 , 79 ). This effect is demonstrated in Fig. 6A and B , where we compare the raw velocity field to the residuals, respectively, with insets that emphasize the importance of correlating the residuals in our directed migratory system. The averaged correlation over time and its statistical comparisons are shown in Fig. 6C and D . Representative correlation curves are shown in Figure S9 (Supplementary Material) , and our approach is described in Methods.

We first analyzed correlation length parallel to the field axis, C x ( r ), as shown in Fig. 6C . For reference, we initially quantified the x -axis velocity correlation length of control tissues over the experimental time period (Fig. 6C , dashed black curve), confirming the expected gradual increase over time ( 63 ). content_tagsBy contrast, during stimulation, we noted a significant increase (nearly 40% by the end of stimulation; Fig. 6C red curve, and center of bar plot) compared to unstimulated control tissues of C x ( r ), followed by a similarly strong reduction relative to control tissues after the stimulus was removed (Fig. 6C , right side of bar plot). However, unlike our other metrics such as the directionality order parameter or mean velocity (e.g. Figs 2 and 3 ), C x ( r ) showed no signs of equilibrating to the control tissue behavior, even 6 h poststimulation—twice the stimulation period (Fig. 6C , plateau trend from 7 h on, and right side of bar plot). As mechanical strain state exhibited different behaviors parallel and perpendicular to the field axis, we also analyzed how correlation lengths changed perpendicular to the field axis, C y ( r ) (Fig. 6D ). Here, we saw a significant decrease in correlation length (approx. 25%, Fig. 6D red curve, and center of bar plot) relative to the control tissue by the end of the stimulation period, followed by similar behavior to what we observed in C x ( r ) poststimulation, where C y ( r ) again remained significantly depressed relative to correlations in an unstimulated control tissue by the 10th h (Fig. 6D , plateau trend from 7 h on, and right side of bar plot).

Electrical stimulation has a lasting effect on spatial correlations. ( A ) Velocity heatmap during electrotaxis showing the global stimulus response of the central 1 × 1 mm region of a representative tissue. Inset shows velocity vector field in a 200 × 100 μm zone. ( B ) Corresponding heatmap of velocity residuals during electrotaxis emphasizing many small, correlated domains. Corresponding inset to (A) showing underlying domains. ( C ), left Temporal dynamics of correlation length for V x : C x ( r ) showing increased correlation during stimulation and decreased correlation poststimulation. (C), right Statistical comparisons of C x ( r ) relative to control (*, P = 0.04) emphasize ∼40% relative increase during stimulation and ∼40% decrease poststimulation. ( D ), left Temporal dynamics of correlation length for V y : C y ( r ). (D), right Statistical comparisons of C x ( r ) relative to control (**, P = 0.004; *, P = 0.04). Note sustained drop in correlation relative to control by 10 h. For all correlations, shading and error bars indicate SD; N stimulated = 9 and N control = 6.

These data emphasize several key details. First, there appears to be a trade-off between C x ( r ) and C y ( r ) during stimulation, where stronger correlations of cell migration fluctuations along the field axis necessitate weaker correlations in the orthogonal axis. More surprising, perhaps, was the clear reduction in correlation length that occurred in all axes poststimulation relative to a control tissue (right sides of bar plots in Fig. 6C and D ). This reduction in correlated domain size again implies that the act of removing the field stimulus appears to trigger a long-term, large-scale loss, or a ‘reset’ of the tissue state relative to an unperturbed tissue—consistent with the resetting of biophysical behaviors shown in Fig. 4 . content_tagsFurther, this emphasizes that our earlier hypothesis was incorrect—the poststimulation speed and directionality relaxation dynamics do not, in fact, reflect a restoration of baseline correlations.

Discussion

This study furthers our understanding of bioelectric manipulation of collective cell migration and general control of group dynamics in 2 key regards. First, we comprehensively assessed the complete step response of a tissue undergoing electrotactic migration, from onset of directed migration to long-term relaxation. Next, we linked these temporal dynamics to specific, supracellular patterns of electrotaxis at the tissue scale. Together, these results revealed both that short-term driven migration can produce longer-term changes to collective cell behaviors (e.g. Figs 3 and 4 ), and that these effects play out differently throughout the tissue (e.g. bulk, leading, trailing, or perpendicular edges) as in Figs 1 , 2 , and 5 . More broadly, our data imply that globally driven directed cell migration (mediated by electrotaxis here) appeared to ‘reset’ key mechanical aspects of the stimulated tissue, both during and long after stimulation, affecting not only collective tissue growth dynamics (Fig. 4 ), but also cell–cell interaction range (Fig. 6 ), apparently mediated by alterations to how mechanical strain coupled across cells (Fig. 5 ).

How migratory programming affects a tissue in time is at least as important as the spatial, or supracellular response, as large variations in response timescales contributed to the overall emergent collective behaviors during and after stimulation. Many tissues, such as epithelia, are tightly cohesive due to cell–cell interactions, which gives rise to both elastic and viscous tissue mechanics at various timescales ( 41 , 80 , 81 ). In our stimulated tissues, rapid responses (< 1 h) include the sharply recoiling edge zones (Fig. 2 ), the short poststimulation speed equilibration (Fig. 3 ), and the quickly propagating and shifting parallel strain waves during and after stimulation (Fig. 5E ). These fast biomechanical responses contrast with the multihour equilibration periods after stimulation, such as the sustained rightward speed drop in recovering edges (Fig. 2C ), the very slow reversion to undirected migration in the tissue bulk (Fig. 3A , C , and E), and the return to coordination through strain waves perpendicular to the field (Fig. 5D – E ). Perhaps most surprising was that coordination through strain waves parallel to the field and cell–cell correlations in the bulk of the tissue (Figs 5D , E and 6C , D ) displayed no sign of recovery during the entire poststimulation period—a regime at least twice as long as the initial stimulation. These are reminiscent of the different timescales in response to mechanical stretch ( 41 , 80 , 81 ), and while electrotaxis does not necessarily stretch the footprint of the tissue as a whole ( 56 ), it does induce displacement and internal deformation by directing tissue flow differently from that observed in an unperturbed tissue. As we have shown that this flow plays out in a tissue with supracellular behavioral zones, it should necessarily cause longer-term plasticity of tissue properties as each zone behaves differently from each other. limitationsFuture studies integrating these results with whole-tissue traction force and cell–cell force analysis would help to clarify this.

Such plasticity and susceptibility to control, along with the apparent mechanical reset (Fig. 4 ), might explain why the tissue maintains such a strong memory of directionality in the tissue bulk (Fig. 3 ) poststimulation, as the underlying mechanical state had been entrained to the command direction. If each and every cell can sense the electrical stimulation, the global command could eliminate the need for tissues to coordinate flow with neighbors. limitationsWhile single-cell mechanical analysis would be needed to better explore this, we saw clues in our correlation length analysis (Fig. 6 ). In the direction of stimulation, fluctuations around the mean velocity undergo an increase in length scale during electrotaxis—implying longer-range cell–cell coordination—which may be necessary to maintain tissue integrity at 2–5× increased migration speeds induced by electrotaxis (Fig. 3B ). Further, the coordination realized during electrotaxis must be qualitatively different than that of standard tissues as it is asymmetric (Fig. 6C vs. Fig. 6D ) and is quickly lost after stimulation is removed. This is consistent with previous work showing that the tissue tension profile aligns perpendicular to the stimulation direction ( 56 ), along with other studies linking the relative strength of E-cadherin mediated cell–cell adhesion with overall susceptibility of a tissue to electrotaxis ( 7 , 25 ).

More generally, that even a simple, universal ‘command’ produced such nuanced behavior highlights not only the importance of taking into account supracellular, or supracollective regional behaviors when controlling groups but, perhaps most importantly, that these behaviors may also differ from each other and from expectations poststimulation. Specifically, our data clearly demonstrated that turning off a stimulus does not simply allow a collective system to return to nominal behavior, meaning that stopping stimulation can be just as transformative to tissue or group dynamics as starting it. These concepts are important to consider in the control of any collective system, spanning theoretical active fluids ( 82 ), large animal groups ( 33 , 83 ), and cells within tissues ( 19 , 38 ). In the longer-term, understanding the relevant timescales and coupling behaviors in driven collective motion will likely help to improve stimulation and control efficacy. Understanding the long-term effectiveness of short-term migration stimulation in a given system can help reduce the need for long-term or more complex stimulation—especially valuable in contexts such as in vivo applications or therapeutics. This is particularly relevant in the growing space of electroceuticals and bioelectronic interfaces being developed to augment large-scale, multicellular processes such as skin healing. Here, there is already striking data of the potential for electric fields to improve re-epithelialization in vitro ( 9 , 84 ) and in vivo ( 23 , 24 , 85–87 ). For instance, we previously demonstrated that closure of the classic in vitro ‘scratch wound’ can be accelerated by at least 2x in vitro using field stimulation ( 9 ), while direct electrical stimulation in vivo of zebrafish tail injuries improved the regeneration process ( 24 ). This is a highly interdisciplinary area of research, where there is an increasing need to investigate how stimulation affects general collective cellular responses. Our work here helps connect the stimulus to both short- and longer-term collective migration responses. Beyond the aforementioned translational applications, the approaches and metrics developed here should also provide a foundation for comparisons of different classes of global commands in tissues (e.g. chemotaxis, mechanical strain, or optogenetics), which may lead to a broader understanding of how to effectively manipulate collective migration more generally. Finally, the distinct spatiotemporal responses we demonstrated and the related biomechanical changes emphasize the likelihood that ‘optimal control’ of collective cell migration will likely require stimuli that are better tuned to specific functional regions within a group.

Materials and methods

Cell culture

MDCK-II wild type and MDCK-II E-cadherin:DsRed ( 68 ) canine kidney epithelial cells were a gift from the Nelson Laboratory at Stanford University and were cultured in customized media consisting of low-glucose (1 g/L) DMEM with phenol red (Gibco; powder stock), 1 g/L sodium bicarbonate (lower than standard DMEM), 1% streptomycin/penicillin (Gibco), and 10% fetal bovine serum (Atlanta Biological). Cells were maintained at 37°C and 5% CO 2 in humidified air.

Micropatterning of epithelial arrays

Standardized arrays of epithelia were produced using silicone tissue stencils built by razor writer as described previously ( 19 , 29 ). A brief summary follows. A total of 250 μm thick sheets of silicone elastomer (Bisco HT6240) were cut into stencil patterns (5 × 5 mm square arrays) using a razor writer (Silhouette Cameo) and transferred into culture vessels. Suspensions of MDCK cells were then seeded into the stencil patterns with the volume and density tuned to produce uniformly dense tissues. A total of 3 ranges of cell densities were analyzed here as follows. Medium density tissues, averaging 2,740 ± 320 (SD) cells/mm 2 , were our baseline samples and intentionally covered a large range. Tissues classified as high density tissues averaged 4,450 ± 90 (SD) cells/mm 2 and exhibited limited movement at the start of the experiment, as was expected ( 60 ). Low density tissues, averaging 2,230 ± 90 (SD) cells/mm 2 , were at the lowest density while still maintaining confluence in the 5 × 5 mm space. In all cases, cells were cultured for 1 h to allow for adhesion and then the dish was flooded with media and maintained for 16 h in an incubator prior to being inserted in our electro-bioreactor as described below.

Electro-bioreactor design

The electro-bioreactor design used here was modified from our SCHEEPDOG platform ( 19 ) for uniaxial stimulation over a large culture area as described previously ( 7 ). Briefly, a custom laser-cut acrylic housing, combined with silicone adhesive layering for a tight seal, was placed around a microarray of tissues in a 10 cm tissue-culture dish to allow media perfusion from north to south and electrical stimulation from west to east. The electro-bioreactor accommodated a 35 × 17 mm cell channel, flanked on both sides by salt bridges of 4% w/v agarose in phosphate-buffered saline (PBS), through which the electric current passed from Ag/AgCl electrodes, each sitting in a reservoir of 2 mL PBS. Titanium probes were inserted into the agarose salt bridges to monitor the voltage across the channel throughout the assay; by connecting these to a USB oscilloscope (Analog Discovery 2, Digilent Inc.), we were able to finely tune the electric current sourced by a Keithley 2450 SourceMeter (Tektronix) with proportional feedback control using a custom MATLAB script to ensure 3 V/cm during stimulation and no current otherwise. For higher resolution images, additional experiments were conducted with a second similarly-modified device with a glass-bottom, coated with 50 μg/mL type-IV collagen (MilliporeSigma). This glass-bottom design was only used for data contributing to Figure S7 (Supplementary Material) and Movie 5, to further accentuate the molecular minutiae during electrotaxis. For a complete summary of the approach and guidelines, see Zajdel et al. ( 19 ) and our repositories ( github.com/cohenlabprinceton ).

In this study, every electrotaxis experiment began by removing the stencil surrounding the tissues 2 h before data collection. After a 1 h unstimulated control period, electrical stimulation was induced at 3 V/cm across the channel for 3 h before the field was removed while imaging continued for the remainder of the experiment. For the assay with chemical junctional disruption, epithelia were otherwise treated as control tissues, but the perfusion media fed into the electro-bioreactor was doped with 4 mM of EGTA (EMD Millipore Corp.) for 30 minutes, from the 3.5 to the 4 h mark, timed to end at the same time as the end of stimulation in our electrotaxed tissue experiments.

Time-lapse imaging and data collection

For MDCK-II wild type experiments, we captured time-lapse microscopy images every 10 minutes with an automated inverted microscope (Zeiss Axio Observer Z1) equipped with a 5x/0.16 phase contrast objective, a Photometrics Prime (Photometrics, Inc.) sCMOS camera, an XY motorized stage, and controlled using Slidebook (3i Intelligent Imaging Innovations). For MDCK-II cells stably expressing E-cadherin:DsRed experiments ( 68 ), we instead used a 20x/0.75 fluorescence objective on an inverted microscope (Nikon Ti2), with a Nikon Qi2 CMOS camera, equipped with an XY motorized stage and controlled using NIS Elements software. In both experimental setups, the microscopes were equipped with a custom-built cage incubator to maintain 37°C, inside which fresh media with continuously bubbled 5% CO2 was perfused with a peristaltic pump (Instech Laboratories) through the electro-bioreactor at a rate of 2.5 mL/h.

Image postprocessing, PIV, and nuclear cell tracking

FIJI ( https://imagej.net/software/fiji ) was used to process all tiled time-lapse images through template-matching ( 88 ), stitching ( 89 ), and masking. Velocity vector fields were calculated from the processed time-lapse images using PIV based on PIVlab ( 90 , 91 ), with a 2-pass iteration of 64 × 64 pixel and 32 × 32 pixel interrogation windows, both with a 50% step size, providing a 16-pixel final step size between vectors. We filtered out vectors that were outside 5 SDs and replaced them with interpolation. All data was then analyzed with MATLAB (Mathworks, 2019a), providing visualizations of the cell movements across the entire tissue.

Cell nuclear and density data were calculated by segmenting images from our in-house Fluorescence Reconstruction Microscopy tool ( 92 ), a convolution neural network trained to produce images of nuclei from our phase contrast images. We also used this output to calculate trajectories of cell nuclei within a 150-μm thick zone from the edges, taking only the center 2 mm near each edge. We linked cell tracks with Linear Motion Tracking, after detecting spots filtered by the Laplacian of the Gaussian, using FIJI’s TrackMate plugin ( 93 ).

Average kymographs

We first created representative kymographs by averaging over the vertical direction across the entire width of each tissue, excluding the outer 1.5 mm of the top and bottom. We then aligned kymographs based on the initial centroid of each tissue and averaged across all tissues for each condition. We did not plot points at the outer boundary that had not yet been reached by more than a single replicate.

Strain wave propagation analysis

PIV analysis requires grayscale images, so we mapped \(\dot{\epsilon }_{xx}\) and \(\dot{\epsilon }_{yy}\) vector fields from all timepoints to 8-bit grayscale image sequences. As such, we mapped the maximum and minimum 0.1% of strain rate values from each vector field to 255 and 0, respectively, interpolating the remaining values to the appropriate integer between 0 and 255. We then performed PIV analysis on the resulting image sequences, using 32 × 32 followed by 16 × 16 pixel windows with 50% overlap. Note that the resulting fields of strain rate propagation will then have resolution 8x less than strain rate. Each pixel of strain propagation then represents a 256 × 256 window of the original phase contrast images.

Characteristic time-scale analysis

Characteristic time scales were obtained by fitting exponential decay curves to mean | V x | and directionality order parameter φ data using a custom MATLAB script. A nonlinear model was fit to the data after the end of stimulation, from the 4th through 10th h, using the model functions \(f(x) = C e^{-\frac{x}{\tau }}\) for mean directionality order parameter and \(f(x) = A + C e^{-\frac{x}{\tau }}\) for mean | V x |. Initial values for the model were taken directly from the data. R 2 was greater than 0.91 for all fits.

Correlation length analysis

Spatial correlation lengths were calculated using a custom MATLAB script. Briefly, 2-D spatial autocorrelation was performed on PIV vector field data utilizing the xcorr2() function, and a radial scan of the autocorrelation matrix was used to obtain a correlation curve. The correlation length was obtained by interpolating the correlation curve and calculating the distance at which the correlation drops below 10% of its maximal value.

Traveling wave analysis

Inward traveling waves of cell mobilization visualized in mean directionality order parameter kymographs of both electrically stimulated and EGTA-treated (chemical junctional disruption) tissues alike were analyzed using a custom MATLAB script. Kymographs were trimmed and binarized before masked using a watershed. Regions of the image were then filled and closed boundaries were located using built-in morphological operations. The traveling waves were then found by finding the first nonzero pixel in each respective half of the image. For speeds, a linear fit was made to the discovered waves, and slopes were reported.

Cellular shape analysis

In order to extract cell boundaries, we first employed our in-house Fluorescence Reconstruction Microscopy tool ( 92 ), as described above, to produce images of nuclei from phase images. Then, FIJI was used for segmenting and processing those images, in particular utilizing ImageJ resident functions ‘Find Maxima’ to resolve a single point for each cell nucleus, and ‘Voronoi’ to construct a Voronoi tessellation of the tissue. This method has been shown to replicate cell membrane tessellation with sufficient accuracy ( 66 , 67 ), as we show with corresponding phase and Voronoi diagram panels in Figure S3 (Supplementary Material) . Individual Voronoi cell perimeter and area were measured using FIJI function ‘Analyze Particles’ and exported to MATLAB for further analysis. Shape index was calculated for each individual cell using the equation (2) \[\begin{eqnarray*} {q} = ( {P}/\sqrt{ {A}}), \end{eqnarray*}\] where P is the perimeter and A is the area of each cell ( 62 , 64 , 65 ). Kymographs of this data ( Figure S6, Supplementary Material ) were then created as described in the ‘average kymographs’ section above. Average cell shape dynamics ( Figure S4, Supplementary Material ) were calculated by averaging the cell shape indices of all cells in a 3 × 3 mm square in the center of stimulated and control tissues, respectively collated within each category of low, medium, and high densities.

Orientation analysis

Orientation analysis was conducted on time-lapse images of single cells, small clusters, and confluent tissues of MDCK-II cells stably expressing E-cadherin:DsRed. Orientation data were produced using ImageJ plug-in ‘OrientationJ’ with a Structure Tensor local window size of 8 pixels (∼2.9 μm) and a cubic spline gradient. Data were stored in HSV images and then imported to MATLAB for further analysis. Data were filtered for background noise before binned and visualized using the polarhistogram() function.

Statistical analysis

Statistical tests were conducted using GraphPad Prism 9.1.2 (GraphPad Software) with an unpaired 2-sided Mann–Whitney nonparametric U test. When comparing a variable between stimulated and control tissues, ensemble values (and not, for example, individual PIV vector data) for each tissue were calculated and then each classification group was compared one against the other.

Supplementary Material

Click here for additional data file.