Cell-cell interactions and fluctuations in the direction of motility promote directed migration of osteoblasts in direct current electrotaxis

Abstract

Under both physiological (development, regeneration) and pathological conditions (cancer metastasis), cells migrate while sensing environmental cues in the form of mechanical, chemical or electrical stimuli. In the case of bone tissue, osteoblast migration is essential in bone regeneration. Although it is known that osteoblasts respond to exogenous electric fields, the underlying mechanism of electrotactic collective movement of human osteoblasts is unclear. Here, we present a computational model that describes the osteoblast cell migration in a direct current electric field as the motion of a collection of active self-propelled particles and takes into account fluctuations in the direction of single-cell migration, finite-range cell-cell interactions, and the interaction of a cell with the external electric field. By comparing this model with in vitro experiments in which human primary osteoblasts are exposed to a direct current electric field of different field strengths, we show that cell-cell interactions and fluctuations in the migration direction promote anode-directed collective migration of osteoblasts.

1 Introduction

The response of the cell to its sensory inputs plays a crucial role in many biological processes such as embryonic development, tissue formation/regeneration, and wound healing. One of the common but crucial reactions of cells is their directed motility, where cells alter their motion in response to external stimuli. Generally, such stimuli are considered to consist of chemical (chemotaxis) or mechanical (adhesion and substrate contact; haptotaxis) mechanisms, as well as of temperature gradients (thermotaxis) or electric fields (electrotaxis) (Simpson et al., 2017), (Lara Rodriguez and Schneider, 2013), (Wu and Lin, 2011), (Zajdel et al., 2020). The latter, electrotaxis, also termed galvanotaxis, is being increasingly studied, in particular in keratinocytes and fibroblasts, since it may provide a promising strategy to foster skin wound healing (Liang et al., 2020), (Cho et al., 2018), (Lin et al., 2017), (Saltukoglu et al., 2015). Electric fields are an attractive approach for manipulating cell behavior because of their broad application across multiple cells and tissue types (Dawson et al., 2020; Gruening et al., 2021). In most cells, the direction of migration seems to be cathodal (e.g., in fibroblasts and mesenchymal stem and corneal epithelial cells (Lin et al., 2017) (Zhao et al., 2002). However, other cell types, such as adenocarcinoma cells, and bone marrow mesenchymal stem cells, show the opposite direction, i.e., anodal (Zhao et al., 2011; Lin et al., 2017). Interestingly, a reversal of directionality was reported for keratinocytes when inhibiting P2Y, a class of extracellular signal sensing surface receptors (Saltukoglu et al., 2015). We recently reported that store-operated calcium channels are pivotal for electrotaxis in human osteoblasts (Rohde et al., 2019), which interestingly migrate to the anode. Thus electrotaxis seems to depend on a variety of factors, such as cell type, environment, possibly age, and ontogenetic stage, all of which should influence cellular signaling pathway equipment.

One of the factors that has not found much consideration so far is that in vivo, electrotactic cell migration involves not only singular but many cells, for example, in a tissue, which collectively responds to either endogenous or exogenous electric fields. Such an electric field-dependent collective cell migration raises the question of how cell interactions, with the electric field on the one hand and the neighbouring cells on the other, shape the final migration vector. In previous modeling studies mainly reaction-diffusion based models were used (Gruler and Nuccitelli, 2000), (Schienbein and Gruler, 1993), (Camley and Rappel, 2017), in some cases cell interaction with the electrical field and chemoattractants were included (Vanegas-Acosta et al., 2012), (Wu and Lin, 2011). These approaches focused on cell migration mainly at the mean-field level and did not resolve the processes at the level of a single cell. Recently, the impact of external electrical field on single-cell motility was quantified using a Bayesian inference approach (Prescott et al., 2021). Cell-cell communication establishes a network that gives rise to many interesting behaviors, such as non-linear collective response, as observed in quorum sensing, a type of bacterial cell-cell communication (Waters and Bassler, 2005; Thurley et al., 2018). Quantitative studies have shown that collective cell migration in epithelial structures is an emergent phenomenon that cannot be explained without considering cell-cell interactions (Barton et al., 2017; Henkes et al., 2020). Self-propelled particle-based models have been widely used to study collective behavior in cell migration in tissues (Szabó et al., 2006; Trepat et al., 2009). While the inclusion of cell-cell interactions in models seems to be natural in the case of high-density tissue culture, where cells adhere to each other and thus exert a pulling force on the neighboring cells, for example, in epithelial wound healing (Brugués et al., 2014), (Cohen et al., 2014), (Zajdel et al., 2020), the rules governing such an interaction in a system of isolated cells, such as in vitro cell culture, remains ambiguous. Three-dimensional multi-cell models based on the finite element method computationally studied the influence of multiple signals, such as the substrate stiffness, chemical gradient, and the electric field on the migratory behavior of cell collective (Mousavi et al., 2013, 2014; Mousavi and Doweidar, 2015; Urdeitx and Doweidar, 2020).

Here, we present a computational model to study the motility and collective migration of osteoblast cells stimulated by direct current (DC) electric field. We compare our model with quantitative experiments and test the hypothesis that cell-cell interactions shape the total vector by re-analyzing data on an individual cell basis from our recent study on osteoblast migration mechanisms in DC electric field (Rohde et al., 2019). Our modeling framework incorporates two-dimensional motion of single cells on a substrate, fluctuations in the direction of cell migration, and key cell interactions. We computationally explore our model to study the influence of two distinct biologically relevant model parameters on electrotactic migration of osteoblasts.

2 Materials and methods

In this study, data on cell migration of human osteoblasts under DC-electrical stimulation were re-analyzed using a previous set of experiments (Rohde et al., 2019). Cell cultivation and stimulation methods are detailed in this paper and given in brief below:

2.1 Cell culture

Human osteoblasts were isolated from femoral heads of patients (n = 14) undergoing a total hip replacement. Patients gave consent and the study was approved by the local ethics committee (permit A 2010–10). Osteoblasts were isolated from cancellous bone as previously described (Lochner et al., 2011). Isolated cells were cultured in Dulbecco’s Modified Eagle Medium (Pan Biotech, Aidenbach, Germany) supplemented with 10% fetal calf serum, 1% amphotericin B, 1% penicillin-strepto-mycin and 1% hepes- buffer under standard cell culture conditions (5% CO₂ and 37°C). Ascorbic acid (50 μg/ml), β-glycerophosphate (10 mM), and dexamethasone (100 nM) (Sigma Aldrich, St. Louis, MO, United States) were added to cell culture medium to maintain osteoblast phenotype. For cell migration experiments cells in passage three were used.

2.2 Direct current electrical stimulation chamber and experimental procedure

To study migration of osteoblasts in electric fields, we used a two-part stimulation chamber described in (Rohde et al., 2019). Before each use, both chamber parts were cleaned with 70% ethanol, washed with a mild detergent and rinsed extensively with distilled water before steam sterilization. Coverslips (24 × 50 mm) for seeding osteoblast cultures were coated with rat tail collagen (Advanced Biomatrix, San Diego, CA, United States) by incubation of 50 μm/ml rat tail collagen diluted in sterile 0.1% acetic acid for 1 h. Coverslips were positioned in a groove in the upper chamber part and edges sealed with silicon paste (Korasilone, Obermeier GmbH, Bad Berleburg, Germany). Upper and lower chamber parts were bolted by 12 screws to ensure tight contact and prevent leakage and chambers were exposed to UV light for sterilization. After this sterilisation treatment, remaining solution was aspired and coverslips were washed twice with phosphate buffered saline (Biochrom, Berlin, Germany) before cell seeding. A total of 2 × 10³ osteoblasts were seeded per chamber and cells were allowed to adhere for 30 min. Afterwards, coverslips were washed twice with medium to remove non-adherent cells. Chambers were then sealed with a top coverglass, and silicon paste and cells accommodated to chamber overnight. For DC-stimulation, silver/silver chloride electrodes were placed into outer reservoirs separated from cell area to avoid electrochemical reactions within the tissue chamber. Current was conducted to the cell chamber using agar bridges (silicon tubes, length 120 mm, inner diameter 5 mm) consisting of 2% agarose (TopVision agarose, ThermoScientific, Waltham, MA, United States) in Ringer’s solution (Braun, Melsungen, Germany). Current was applied to the electrodes for 7 h via crocodile clamps using a DC power supply (Standard Power Pack P25, Biometra, Göttingen, Germany). To maintain constant stimulation, voltage was measured directly at the borders of the cell area (electrode distance 24 mm) using a multimeter (Voltcraft VC220, Conrad Electronic, Wollerau, Switzerland) and adjusted during the experiments. Each of the experiments was conducted with one cell culture being divided to obtain a sham stimulation group as control, and a DC-stimulation group for the respective field strength used. Electric field strengths were 160, 300, 360, 426 and 436 V/m.

2.3 Migration analysis

For the analysis, all cells from the sham groups were pooled as one control. Thus, a total of n = 177 (sham), 34 (160 V/m), 35 (300 V/m), 26 (360 V/m), 43 (426 V/m) and 33 (436 V/m) cells were analysed. For this, photographs were taken at eight fields of view evenly distributed over the cell area at beginning (Figure 1A) and end time points (Figure 1B) with a Leica DMI 6000 and LAS X software during the 7-hour (h) stimulation, or sham stimulation, procedure. The pairs of photographs were then aligned manually, and merged, taking external markers as reference points (Figure 1C). To quantify migration within the electric field, segmentation of the cell shape, including cell extensions, was performed manually using Image J software (NIH) for each cell that could be identified in both the time points (see yellow coastlines in Figure 1C), i.e., 0 h and 7 h after DC stimulation. The segmentation procedure yielded, among other features, the coordinates of the cell centroid. Using the coordinates of the cell centroid at these two time points, i.e., 0 h and 7 h, the distance and direction of migration was calculated for each cell. The migration distance was defined as , where X1, Y1 and X2, Y2 represent the coordinates of the cell centroid at 0 h and 7 h after DC stimulation, respectively Figure 1C. The migration angle was defined as . Using the distance and direction of cell migration we obtained a migration plot for each case of electrical stimulation, which could be depicted in a polar coordinate system, as shown in. The anode in the polar plots of DC stimulated experiments is located at 180° angle. For better comparison of all experiments, we binned the migration angles into 36 sectors of 10° each, and classified migration speeds in a scoring system. Thus, the migration angle was calculated starting from the original cell position, and angles were assigned to the 36 sectors, where sector 10–18 (90–180°) and 18–26 (180–270°) represent anode-directed migration, while sectors 1–9 and 27–36 (0–90 and 270–360°) represent cathode-directed cell migration. To construct polar plots (Figures 1D,E) illustrating both migration direction and velocity, the migration speed of single cells was colour coded from 0 to 18 μm/h in nine groups of 3 μm/h bins. The relative sector lengths denote the percentage of cells migrating at a certain speed range.

FIGURE 1

3 Results

3.1 In vitro direct current stimulation of human osteoblasts

In the experimental part of this study, we exposed human osteoblasts to DC electric fields for 7 h at different stimulation strengths and matched each of these experiments with a sham-stimulated, control group treated identically, save the DC stimulation. For the analysis of the migration behaviour, we selected adherent cells in the stimulation chambers which could be identified clearly at starting and end time points of the experiment, and did not form clusters precluding the outlining of their boundaries (Figures 1A–C). Using photographs of several fields of vision in each chamber, 1–4 cells could be traced in this way per field of vision position, totalling n = 177 cells (sham stimulation), as well as n = 34 cells (at 160 V/m), n = 35 cells (at 300 V/m), n = 26 (at 360 V/m), n = 43 (at 426 V/m) and n = 33 (at 436 V/m). The center position coordinates were determined for each cell tracked between the two time points, and this was used to calculate the effective magnitude and direction of the cell velocity.

As one can notice in the original photographs of one typical cell from the experiment using 436 V/m stimulation, the cells move, in this case anodally, and at the same time change their shape within the 7 h stimulation (Figures 1A–C). While we did not analyse shape changes any further in this study, we took them into consideration for determining the centroids of the cells (colored dots in Figure 1C), which we used as markers to determine the net movement (blue dashed arrow in Figure 1C).

Comparing cell migration velocities (plotted as sectors of polar plots) without stimulation (Figure 1D), to those with weak (160 V/m; Figure 1E) or strong stimulation (436 V/m; Figure 1F), one can appreciate that the directionality of migration shifts with increasing field strength from random, covering all sectors of the plot (Figure 1D) to exclusively anodal, covering only the anodal sectors (Figure 1F). At the same time, also the speed of the cells appears to shift from lower speeds with a maximum of 8–10 μm/h (green hues in Figure 1D corresponding to of the cells) to a maximum of 16–18 μm/h (red sectors in Figure 1F corresponding to of the cells). To address the question of a possible correlation of speed and field strength, we quantified the cell migration speed of all cells in all experiments under different stimulation strengths. As shown in Figure 1G, migration speed under DC-stimulation is significantly different from sham stimulation conditions only at 300 V/m and 436 V/m (p 0.001, ANOVA on ranks with Dunn’s all-pairwise comparisons). We further performed a linear regression on the cumulative data for all the field strengths (blue line in Figure 1G) which showed a weak correlation between field strength and migration speed, with a regression coefficient of R = 0.516.

3.2 Modeling electrotactic migration of osteoblast cells

Our computational model simulates the in vitro motility behaviour of individual cells that are subject to external DC electric field, (Figure 2A). The main components of our model are (i) the ability of the cells to interact with the other cells, (Figures 2B,C), and, (ii) the ability of the cells to interact with the external electric field, (Figure 2D). Cell-cell interaction involves two types of forces: short-range repulsive forces and the alignment of the direction of motion with the cells’ local neighbors, (Figures 2B,C). Through soft-core repulsion, the force at short distances ensures that cells do not overlap. We also include in our model the influence of each cell’s local neighbors on the direction of its migration. Such cell-to-cell interactions certainly play a role in high-density tissue culture via cell-cell contacts. However, since mechanical or signaling cues are at least conceivable also in 2D cell cultures without direct cell contacts, we introduce this factor in the model to study the possible role of such interactions in our experiments. Finally, we also consider the interaction of cells with the applied DC electric field, (Figure 2D). We also include angular noise that describes the uncertainty that a cell experiences in directional sensing during migration, (Figure 2E).

FIGURE 2

3.2.1 The cell

We describe each cell as a circular disk of radius R = 1 exhibiting overdamped dynamics in two dimensions (2D) with an active speed of v₀. The state of each cell i is characterized at time t by its position , described through the coordinates , and its migration velocity , where, v₀ is the cell migration speed and is the unit vector representing the direction of cell migration, with θ being the angle that the cell makes with the horizontal axis of the laboratory frame. The direction θ that a cell takes at the next time step depends not only on its direction of motion in the immediately preceding time step, but also on the forces acting on the cell. The total force acting on the cell i results from cell-cell interactions and cell interaction with the applied DC electric field. These forces are discussed in more detail in the following sections.

3.2.2 Cell-cell interactions

We consider two types of cell-cell interactions in our model, one due to finite-volume exclusion and the other resulting in the alignment of a cell’s migration direction with its neighbors. Each cell is assumed to occupy a finite area in the cell culture medium in which it is placed. In order to ensure that the cells do not overlap, we include, at each time t, a repulsive force K_ij that is proportional to the degree of overlap between any two cells and is given by,where, is the euclidean distance between the center of two cells i and j at time t and is calculated as . is the unit vector pointing from the center of cell i to the center of cell j. The total repulsive force acting on cell i at time t, denoted by is the sum of all the pairwise repulsive interactions between cell i and its nearest overlapping neighbors j,

3.2.3 Cell-electric field interaction

The electric field in this model is defined on a regular square lattice underlying the domain in which the cells are migrating. Each grid point k whose position q_k is specified by its coordinates (m_k, n_k). The electric field at the grid point q_k located at (m_k, n_k) is characterized by the vector E_k = E₀e_k, where, e_k = (cos Θ_k, sin Θ_k) is the unit vector representing the orientation of the electric field, and, E₀ is the magnitude of the electric field which corresponds to the field strength of the electrical stimulation in the experiments. Θ_k is the angle that the electric field vector at (m_k, n_k) makes with the horizontal x − axis. Due to cell spreading on the substrate, the electric field sensors on the cell membrane relocate and are distributed over a finite region in the motility domain of each individual cell. Therefore, we assume that the contribution to the total electrical force acting on a cell comes from the electric field at all the grid points within a region that the cell effectively covers due to cell spreading. We denote the extent of this region by radius r_e taken from the center of each cell i. Thus, the electric field that the cell i experiences is the sum total of the electric field at all the grid points that lie in the region within the radius r_e from the center of cell i and is given by,Assuming that the cell motion is described by overdamped dynamics, the position of each cell i is updated in the next time step, after calculating the resultant displacement and angle θ due to all the interactions, including the directional alignment with the cells’ neighbour, by the following scheme:where, ν is a friction factor that is associated with the cell-substrate interaction, μ is the effective electrical mobility of the cell. The time step is equal to 1.

Finally, in addition to the various forces acting on a cell, as described above, the cell in this model also undergoes directional alignment with its proximal neighbours. The directional alignment of a cell with its neighbours influences the direction of the cell’s motility and is given by,The directional alignment in this model is only hampered by an angular white noise uniformly distributed in with and and whose strength is given by η. arg [c] in the first term in Eq. 5 refers to the angle associated to the vector c if this is expressed in polar coordinates, and the sum is taken over all particles j within distance r_a of cell i (including i itself).

3.2.4 Simulation details

We simulate the motility behavior of N = 35 cells, as in our experiments, there are approximately 30–40 cells in a single field of view. Cells are initially randomly distributed in a circular region within the spatial domain representing the stimulation chamber. Osteoblast cells are roughly 100 μm in diameter, considering all cells extensions, and, we use this to define the cell radius R, which is one length unit in our simulations. The cell radius R is assumed to be the basic length scale in these simulations. Time steps are separated by Δt which is set to 1. The parameters of the model and their values used in these simulations are listed in Table 1. The active speed of cells is 0.1R per time step.

The time parameters in our simulations are scaled such that the speed of the unstimulated cells corresponds to the average cell speed in the experimental sham case, which is ∼2 − 3μm/h. The magnitude of the electric field E₀ = 0.014 corresponds to 160 V/m in experiments and is chosen such that the cell speed, in this case, is approximately equal to the experimental value, i.e., . The magnitude of the electric field corresponding to 436 V/m in experiments is obtained by simply multiplying 0.014 by the factor 2.75, as 436/160 = 2.75. At the start of each simulation, we specify the initial positions x_i (t = 0), y_i (t = 0), initial speed v_i (t = 0) and the velocity direction θ_i (t = 0) of each cell i. Initially, the cells are placed at random locations within the domain. The initial velocity direction is also distributed randomly in the range [0,2π]. At each time step for each cell, we identify cells that are less than a distance of 2R apart. We developed a distance sorting algorithm to efficiently provide a list of nearest neighbors in contact with each cell for use in the simulations. From this we calculate the total repulsive force due to volume exclusion acting on each cell from all of its neighbouring cells, as given by Eqs. 1, 2. Using the same distance sorting algorithm, we also determine all the grid points of the underlying grid, on which the electric field is defined, that lie within a radius r_e of each cell and calculate the net electric field, denoted by D_i as given by Eq. 3. This constitutes the net force due to the electric field D_i acting on each cell i. Experiments show that the osteoblast cell migration is anode-directed. We incorporate anode-directed motility of individual cells into our model by assigning a polarity p, which can assume the value + 1 (cathode-directed) or − 1 (anode-directed), to the net electrical force experienced by the cell as, . In addition, we also determine for each cell all its neighbouring cells that are located within the radius r_a, and calculate the net direction of velocity. The direction of velocity of each cell is then updated by its net direction, to which a weak noise η = 0.05 is added, as given by Eq. 5. Finally, the positions of each cell are then updated using the Eq. 4.

3.2.5 Migratory behavior of osteoblasts in direct current electrical field

To study the influence of an externally applied DC electrical field on human osteoblast cell migratory behavior, we simulated N = 35 migrating cells with and without DC electrical stimulation for 140 time steps. We developed the computational model in (version 2021b) MATLAB (2021), wherein we described the migration dynamics of each cell by the equations of motion Eqs. 4, 5. We use periodic boundary conditions.

We have verified, through multiple runs of the simulations, that our results are qualitatively invariant of random initial conditions and stochastic angular fluctuations in the simulations [(see Supplementary Material)]. Figure 3 left column shows the position of all the cells at the final time step (gray circles with arrows) as well as their trajectories at each time step, in the case of no electrical stimulation, Figure 3A, and in the case of DC stimulation with the magnitude of the electric field being 0.014 (160 V/m) and 0.038 (436 V/m), Figures 3B,C, respectively. The magnitude 0.038 of the electric field in simulations corresponds to the maximum field strength of the electrical stimulation in experiments, i.e., 436 V/m, Figure 1G. The velocity of cell migration, calculated from the initial and the final time step, is shown in Figures 3D–F (right column) as polar plots for the case without electrical stimulation, Figure 3D, and with electrical stimulation of different magnitudes of electric field, i.e., 0.014 and 0.038, Figures 3E,F, respectively. Each polar plot shown in Figures 3D–F, is the cumulate of ten separate runs of the simulation. The initial velocity of each cell and the noise in the direction of cell velocity at each time step are random, and this renders robustness to the results of polar plot distributions.

FIGURE 3

In the absence of electrical stimulation, the cells move, as expected, in all directions, Figure 3A. Trajectories of individual cells show that, over time, all the cells collectively explore the space isotropically, Figures 3A,B feature that is also reflected in the polar plots, 3d. The mean cell speed in this case is m/h. However, when a DC electrical field of magnitude 0.014, which corresponds to 160 V/m, is applied, the cells start exhibiting a directional migration towards the anode 3b. The individual cell trajectory plot shows that although the final position of the majority of the cells is towards the anode, few cells still migrate towards the cathode, albeit at much shorter distances than the anodally migrated cells, Figure 3B. The polar plot, showing the effective cell migration velocity, clearly presents the modulation of the direction of migration by an external field, Figure 3E. Following the trajectories of individual cells also reveals that cell migration is not instantaneously switched in the direction of the anode. Cells respond to the applied electrical field by gradually changing the direction of their migration. Initially, most of the cells move randomly. However, at later times, under the influence of the electric field they gradually turn towards the anode. This delayed response in eventual anode-directed motility of cells is because the electrotactic migration speed μD, where μ is the cell mobility, is much weaker than the active cell migration speed v. Cells display a highly directional motion towards the anode with increasing strength of the electric field E₀, Figure 3C. Cell migration, in this case, shows a much faster re-orientation and more persistent motion towards anode 3c. Figure 3F shows that not only the direction of the motion but the effective speed of cell migration is also influenced by increasing the strength of the electric field. The maximum cell speed in this case even reaches up to 10–12 μm/h, Figure 3F.

To better quantify the changes in the collective cell migratory behavior, we calculate, in both experiments and simulations, the directionality order parameter Φ, which reflects how well cell movements have aligned with the electric field and directed towards the anode and is given by,where N is the total number of cells, and the sum is over the cosine of the migration direction of individual cells θ_i. Since the electric field is aligned along the x − axis, we use the cosine component of the cell velocity to evaluate how well the cells are aligned with the external electric field. Φ can vary between 1 (perfectly aligned towards the cathode) and −1 (perfectly aligned towards the anode), and Φ ≃ 0 corresponds to both random isotropic cell movement and movement perfectly orthogonal to the direction of the applied electric field. Our results show that for the listed choice of parameters, as given in Table 1, the directionality order parameter Φ obtained from the model simulation matches very closely with the experiments (Figure 3G). Our results shown in Figure 3 reproduce the following experimental observations: (i) in the absence of electrical stimulation, which corresponds to the experimental sham case, directional migration of cells is not observed, i.e., cells do not move in a preferred direction and are not aligned with the external electric field, (ii) alignment of the direction of cell migration with the external electric field and directed towards anode depends on the strength of the applied electrical field, (iii) the effective cell velocity increases with the electric field strength.

3.2.6 Influence of radius of alignment and cell number density on cell migration

In our model, the collective behavior results from the directional alignment of individual cells with their neighbors. This is controlled by the model parameter r_a, the distance over which the cell aligns its direction of migration with its neighbors, and η, which is the strength of the fluctuation in the direction of the migration of an individual cell. In the simulation results discussed in the preceding section, Figure 3, we set r_a = 2R, which implies the direction alignment of velocity occurs only when cells touch each other. We wondered how the distance over which a cell can align with its neighbors and the noise in the directional sensing influences the collective behavior of osteoblast cell migration. Therefore, we performed a parameter sweep study of the model parameter r_a, which determines the distance over which a cell aligns its direction of motility with its neighbours, with fixed electrical stimulation of strengths 0.014 and 0.038, which correspond to 160 V/m and 436 V/m in experiments, respectively. Our results show that even in the case of weak electrical stimulation , which corresponds to 160 V/m in the experiments, with increasing r_a the cells move in a more directed manner towards the anode, i.e., Φ approaches the value of − 1, as shown in (Figure 4A). This trend, in the dependence of Φ on r_a, is also observed when the electric field strength was varied, from (160 V/m) to (436 V/m), keeping all the other model parameters same in the two cases (Figure 4B). We also found that by only changing the electric field strength , the model best fits the directionality obtained from experiments for r_a = 2 (Figures 4A,B, where the directionality value obtained from experiments is shown as a green dot). The values of all the other parameters except r_a and are as mentioned in Table 1. This parameter sweep-based analysis, in which only the electric field strength is varied and all the other parameters of the model are kept the same in the two cases, suggests that the applied electric field has the strongest influence on the directionality of osteoblast cell migration, and that the cell-alignment interaction in individual osteoblast cell migration is mainly due to direct contact between the cells. The influence of radius of alignment on the directionality of migration further raised our interest in exploring its possible duality with the cell number density and its influence on the directionality of osteoblast cell migration. The frequency of cell-cell interactions, or contacts, is expected to be higher with increasing cell number density. Using the model, and by performing a parameter sweep analysis, we studied the influence of cell number density on the directionality of cell migration. We considered four different cell densities of 0.02, 0.03, 0.04 and 0.06, each of which correspond to a total number of 25, 35, 45 and 65 cells, respectively. Our results show that the alignment of osteoblast cell migration with the external electric field, and its directedness towards anode increases with cell number density (Figure 4C). For the same values of all the model parameters and an electric field strength (which corresponds to 436 V/m), our simulation results revealed that the directionality Φ approaches a value of − 1 in the case of cell density of 0.06 (Figure 4C).

FIGURE 4

3.2.7 Influence of noise on cell migration

Osteoblast cells do not migrate along a straight trajectory over long time intervals. Using our model, we studied the influence of fluctuations in the direction of migration on both the directionality and the effective velocity of electrotactic migration of osteoblasts. Cell movement showed higher directedness with increasing noise strength η (Figure 4D). Although cells move in a more directed manner towards the anode, the effective cell migration velocity, on the contrary, decreases with increasing noise strength (Figures 4D–F). The plot of cell trajectories shows that although cells move in an increasingly directed manner towards the anode, they, traverse shorter distances as angular fluctuation increases. These results suggest that the parameters, i.e., r_a and η, can significantly alter the dynamics of cell migration and give rise to collective behavior in the electrotactic motion of osteoblast cells.

4 Discussion

The migration of osteoblasts, which plays a key role in bone regeneration, can be modulated by external electrical stimulation (Ferrier et al., 1986). This offers an attractive approach toward building electrically active implants for effective tissue regeneration (Brighton et al., 1985; Kaivosoja et al., 2015; Hiemer et al., 2016). In the present paper, we presented a computational model to study (i) the migratory behavior of osteoblasts and, (ii) the consequences of the application of external electrical fields on their migration. The model was used to study the collective behavior of many cells in in vitro experiments where primary human osteoblasts placed in an electrotaxis chamber were stimulated by a DC electric field. For this purpose, we re-analyzed the galvanotactic migration of human osteoblasts exposed to DC-electric field stimulation at different field strengths from a previous study published in (Rohde et al., 2019), now using single-cell rather than clustered data. As observed in our previous paper (Rohde et al., 2019), we confirmed that field exposition leads to migratory directionality towards the anode, and elucidate that the migratory speed distribution ranges from 2 to 18 μm/h, with significantly higher speeds of migration than unstimulated cells at DC-field strengths of 300 and 436 V/m. At this point, it remains elusive why there is a significant impact on migration speed at field strengths of 300 and 436 V/m, and not at other stimulation strengths, beyond the general observation that non-linearities do often emerge in biological systems. This point certainly merits further molecular exploration. Using this single-cell analysis approach, beyond our initial findings in the cited paper using pooled data, i.e., stimulated vs. unstimulated only, we show quantitatively that the directionality of cell migration is thus actually influenced by the field strength, with random migration without stimulation, ∼ 65% anodal migration at low (160 V/m) and exclusively anodal migration at highest field strength (436 V/m). Our detailed cell-by-cell analysis also shows that, although the directionality of cell migration clearly correlates with the strength of the applied electric field, there is only a weak correlation between migratory speed and electric field strength, a correlation which could not be seen in the pooled analysis of our previous paper.

To explain these experimental observations, we modeled each cell as an active agent whose movement is influenced by its interactions with other cells, external electric field, and stochastic switching in the direction of migration. The model considers the force experienced by the cell due to the applied DC electric field. We also considered two types of intercellular interactions: in addition to the nearest-neighbor interaction that ensures finite-volume exclusion by penalizing cell overlaps, cells also interact with other cells via a velocity alignment mechanism. Although specific molecular mechanisms underlying these interactions remain unclear, two important questions can be addressed by the current simulation study: (i) Does directionality also depend on interaction among neighbouring migrating cells, and if so, how large is this interaction radius, (ii) Do directionality and migration speed depend on the accuracy of the putative cellular field sensing mechanism, i.e., in which way does a noise factor influence migration directionality and migration speed?

Our results show that the motility behavior of cells is influenced by the distance over which the cell aligns with its neighbors, stochastic switching in the direction of migration, and the strength of the applied electric field. The simulations in the present paper closely match the experimentally observed weak correlation between migration speed and the applied electric field and are more realistic than previously published ones (Vanegas-Acosta et al., 2012), which predicted speed ranges from 1.8 to 4.0 μm/s, i.e., nearly tenfold the maximum observed by us. As discussed previously, migration at such high speeds probably finds its limitations in adhesive forces acting on the cells on the one hand and rate-limiting factors such as actin conformational change being limited by temperature and Ca²⁺ dynamics (Sich et al., 2010; Jacobs et al., 2011). We performed a quantitative comparison of the directionality order parameter obtained from simulations with experiments as shown in Figure 3G, where the directionality Φ obtained from experiments practically overlaps with the values obtained from simulations. As the simulation results show, varying r_a from 0 (i.e., the case with no inter-cellular interactions) to 8 (i.e. the case with inter-cellular interactions between two cells extending to distances of four cell diameters), the directionality obtained from our experiments best matches the simulation results in which r_a has the value equal to 2. These results suggest that the interactions between cells only in direct contact likely lead to parallel anodal movement. The mechanism of this interaction could be speculated to rely on e.g., osteoblast binding via cadherin, an interaction known to be important for morphogenesis of osteoblasts and subsequent modulation of actin function (Stains and Civitelli, 2005; Stains et al., 2019). Long-distance effects, mediated by e.g. molecules secreted from the cells, tension changes within the collagen coating or distortion of the electric field by the neighboring cell are, in turn, unlikely to be important for osteoblasts. Using our model we also performed a quantitative parameter sweep study to explore the influence of cell number density on the migration directionality parameter Φ. The computational results show that higher cell number density gives rise to higher overall directedness in the cell migration. At higher cell number density one would expect higher number of cell-cell interactions leading to directional alignment of migration, as observed in our simulations.

Our results show that stochastic orientational switching can significantly alter cellular electrotactic motility behavior. In this case, a perfectly directed motion towards the anode is achieved for very high fluctuation strengths, which appears to be counter-intuitive since one would expect that the accuracy of directional movement aligned with the electric field decreases for higher angular fluctuations. Varying η in our simulations from 0 to 0.8, the directionality of ∼ −0.45 in our experiments is in line only with a very low degree of noise (around 0.05, which corresponds to fluctuations of ° in the direction of cell migration), but not commensurate with values of 0.25. The experimental migration speed found to be in the range of 2–12 μm/h would also cover the simulated value of ∼ 4.75 μm/h at η = 0.05. It is, however, conceivable that other cell types do show more influence of noise (arguably reflecting e.g., less mechanical interactions with the substrate, varying cell shape influences, or different field sensing or signaling mechanisms). What remains to be explained is the seemingly paradoxical result that higher fluctuation levels should lead to higher accuracy in directionality. Our hypothesis would be that higher fluctuation raises the probability of cell-to-cell interactions, which will lead to common field alignment. If this hypothesis holds, such movement will lead to cell alignment of migration direction with the field with higher accuracy but lower speed due to frequent corrective movements. Although experiments are needed to validate this hypothesis, it is interesting to note that at the highest stimulation strength of 436 V/m, those cells which are best aligned to the field and directed towards the anode do not belong to the fastest subset of cells (which are, indeed, 10°–30° off the “ideal” orientation; see Figure 1F.

5 Conclusion

The computational model presented here provides a framework for studying in vitro DC electrotactic migration of osteoblast cells in two dimensions and elucidating the rules and the role of individual cell interactions with other cells and with their physical environment. This model is also relevant to study the influence of additional factors on cell migration, such as the cell density and other modes of electrical stimulation, e.g. alternating current stimulation. The model we present here allows for easy integration of additional details, as more data becomes available. Our approach could serve as a tool to not only test existing hypotheses of electrotactic cell migration but also predict migratory behaviour under perturbation conditions, and thus bridge the gap between single cell and collective response in a more effective manner.

References

• 1

  BartonD. L.HenkesS.WeijerC. J.SknepnekR. (2017). Active Vertex Model for cell-resolution description of epithelial tissue mechanics. PLoS Comput. Biol.13, e1005569. 10.1371/journal.pcbi.1005569

  • CrossRef
  • Google Scholar

• 2

  BrightonC. T.HozackW. J.BragerM. D.WindsorR. E.PollackS. R.VreslovicE. J.et al (1985). Fracture healing in the rabbit fibula when subjected to various capacitively coupled electrical fields. J. Orthop. Res.3, 331–340. 10.1002/jor.1100030310

  • CrossRef
  • Google Scholar

• 3

  BruguésA.AnonE.ConteV.VeldhuisJ. H.GuptaM.ColombelliJ.et al (2014). Forces driving epithelial wound healing. Nat. Phys.10, 683–690. 10.1038/NPHYS3040

  • CrossRef
  • Google Scholar

• 4

  CamleyB. A.RappelW.-J. (2017). Physical models of collective cell motility: From cell to tissue. J. Phys. D. Appl. Phys.50, 113002. 10.1088/1361-6463/aa56fe

  • CrossRef
  • Google Scholar

• 5

  ChoY.SonM.JeongH.ShinJ. H. (2018). Electric field-induced migration and intercellular stress alignment in a collective epithelial monolayer. Mol. Biol. Cell29, 2292–2302. 10.1091/mbc.E18-01-0077

  • CrossRef
  • Google Scholar

• 6

  CohenD. J.James NelsonW.MaharbizM. M. (2014). Galvanotactic control of collective cell migration in epithelial monolayers. Nat. Mat.13, 409–417. 10.1038/nmat3891

  • CrossRef
  • Google Scholar

• 7

  DawsonJ.LeeP.van RienenU.AppaliR. (2020). A general theoretical framework to study the influence of electrical fields on mesenchymal stem cells. Front. Bioeng. Biotechnol.8, 557447. 10.3389/fbioe.2020.557447

  • CrossRef
  • Google Scholar

• 8

  FerrierJ.RossS. M.KanehisaJ.AubinJ. E. (1986). Osteoclasts and osteoblasts migrate in opposite directions in response to a constant electrical field. J. Cell. Physiol.129, 283–288. 10.1002/jcp.1041290303

  • CrossRef
  • Google Scholar

• 9

  GrueningM.DawsonJ.VoelknerC.NeuberS.FrickeK.van RienenU.et al (2021). Automatic actin filament quantification and cell shape modeling of osteoblasts on charged ti surfaces. Appl. Sci.11, 5689. 10.3390/app11125689

  • CrossRef
  • Google Scholar

• 10

  GrulerH.NuccitelliR. (2000). The Galvanotaxis response mechanism of keratinocytes can be modeled as a proportional controller. Cell biochem. Biophys.33, 33–51. 10.1385/cbb:33:1:33

  • CrossRef
  • Google Scholar

• 11

  HenkesS.KostanjevecK.CollinsonJ. M.SknepnekR.BertinE. (2020). Dense active matter model of motion patterns in confluent cell monolayers. Nat. Commun.11, 1405. 10.1038/s41467-020-15164-5

  • CrossRef
  • Google Scholar

• 12

  HiemerB.ZiebartJ.Jonitz-HeinckeA.GrunertP. C.SuY.HansmannD.et al (2016). Magnetically induced electrostimulation of human osteoblasts results in enhanced cell viability and osteogenic differentiation. Int. J. Mol. Med.38, 57–64. 10.3892/ijmm.2016.2590

  • CrossRef
  • Google Scholar

• 13

  JacobsD. J.TrivediD.DavidC.YengoC. M. (2011). Kinetics and thermodynamics of the rate-limiting conformational change in the actomyosin V mechanochemical cycle. J. Mol. Biol.407, 716–730. 10.1016/j.jmb.2011.02.001

  • CrossRef
  • Google Scholar

• 14

  KaivosojaE.SariolaV.ChenY.YtK. (2015). The effect of pulsed electromagnetic fields and dehydroepiandrosterone on viability and osteo-induction of human mesenchymal stem cells. J. Tissue Eng. Regen. Med.9, 31–40. 10.1002/term.1612

  • CrossRef
  • Google Scholar

• 15

  Lara RodriguezL.SchneiderI. (2013). Directed cell migration in multi-cue environments. Integr. Biol.5, 1306–1323. 10.1039/c3ib40137e

  • CrossRef
  • Google Scholar

• 16

  LiangY.TianH.LiuJ.LvY. L.WangY.ZhangJ. P.et al (2020). Application of stable continuous external electric field promotes wound healing in pig wound model. Bioelectrochemistry135, 107578. 10.1016/j.bioelechem.2020.107578

  • CrossRef
  • Google Scholar

• 17

  LinB. J.TsaoS. H.ChenA.HuS. K.ChaoL.ChaoP. H. G. (2017). Lipid rafts sense and direct electric field-induced migration. Proc. Natl. Acad. Sci. U. S. A.114, 8568–8573. 10.1073/pnas.1702526114

  • CrossRef
  • Google Scholar

• 18

  LochnerK.FritscheA.JonitzA.HansmannD.MuellerP.Mueller-HilkeB.et al (2011). The potential role of human osteoblasts for periprosthetic osteolysis following exposure to wear particles. Int. J. Mol. Med.28, 1055–1063. 10.3892/ijmm.2011.778

  • CrossRef
  • Google Scholar

• 19

  MATLAB (2021). Matlab. Natick, Massachusetts: The MathWorks Inc.

  • Google Scholar

• 20

  MousaviS.DoweidarM.DoblareM. (2013). 3D computational modelling of cell migration: A mechano-chemo-thermo-electrotaxis approach. J. Theor. Biol.329, 64–73. 10.1016/J.JTBI.2013.03.021

  • CrossRef
  • Google Scholar

• 21

  MousaviS.DoweidarM.DoblareM. (2014). Computational modelling of multi-cell migration in a multi-signalling substrate. Phys. Biol.11, 026002. 10.1088/1478-3975/11/2/026002

  • CrossRef
  • Google Scholar

• 22

  MousaviS. J.DoweidarM. H. (2015). Three-dimensional numerical model of cell morphology during migration in multi-signaling substrates. PLoS ONE30, e0122094. 10.1371/journal.pone.0122094

  • CrossRef
  • Google Scholar

• 23

  PrescottT. P.ZhuK.ZhaoM.BakerR. E. (2021). Quantifying the impact of electric fields on single-cell motility. Biophysical J.120, 3363–3373. 10.1016/j.bpj.2021.06.034

  • CrossRef
  • Google Scholar

• 24

  RohdeM.ZiebartJ.KirschsteinT.SellmannT.PorathK.KühlF.et al (2019). Human osteoblast migration in DC electrical fields depends on store operated Ca2+-release and is correlated to upregulation of stretch-activated TRPM7 channels. Front. Bioeng. Biotechnol.7, 422. 10.3389/fbioe.2019.00422

  • CrossRef
  • Google Scholar

• 25

  SaltukogluD.GrünewaldJ.StrohmeyerN.BenschR.UlbrichM. H.RonnebergerO.et al (2015). Spontaneous and electric feld-controlled front-rear polarization of human keratinocytes. Mol. Biol. Cell. 10.1091/mbc.E14-12-1580

  • CrossRef
  • Google Scholar

• 26

  SchienbeinM.GrulerH. (1993). Langevin equation Fokker-Planck equation and cell migration. Bull. Math. Biol.55, 585–608. 10.1007/BF02460652

  • CrossRef
  • Google Scholar

• 27

  SichN. M.O’DonnellT. J.CoulterS. A.JohnO. A.CarterM. S.CremoC. R.et al (2010). Effects of actin-myosin kinetics on the calcium sensitivity of regulated thin filaments. J. Biol. Chem.285, 39150–39159. 10.1074/jbc.m110.142232

  • CrossRef
  • Google Scholar

• 28

  SimpsonM. J.LoK. Y.YsS. (2017). Quantifying the roles of random motility and directed motility using advection-diffusion theory for a 3T3 fibroblast cell migration assay stimulated with an electric field. BMC Syst. Biol.11, 39. 10.1186/s12918-017-0413-5

  • CrossRef
  • Google Scholar

• 29

  StainsJ. P.CivitelliR. (2005). Cell-to-cell interactions in bone.

  • Google Scholar

• 30

  StainsJ. P.FontanaF.CivitelliR. (2019). “Intercellular junctions and cell-cell communication in the skeletal system,” Principles of bone biology. 2020, 432–442. 10.1016/B978-0-12-814841-9.00018-X

  • CrossRef
  • Google Scholar

• 31

  SzabóB.SzöllösiG. J.GönciB.JurányiZ.SelmecziD.VicsekT. (2006). Phase transition in the collective migration of tissue cells: Experiment and model. Phys. Rev. E74, 061908. 10.1103/PhysRevE.74.061908

  • CrossRef
  • Google Scholar

• 32

  ThurleyK.WuL. F.AltschulerS. J. (2018). Modeling cell-to-cell communication networks using response-time distributions. Cell Syst.6, 355–367.e5. 10.1016/j.cels.2018.01.016

  • CrossRef
  • Google Scholar

• 33

  TrepatX.WassermanM. R.AngeliniT. E.MilletE.WeitzD. A.ButlerJ. P.et al (2009). Physical forces during collective cell migration. Nat. Phys.5, 426–430. 10.1038/nphys1269

  • CrossRef
  • Google Scholar

• 34

  UrdeitxP.DoweidarM. H. (2020). A computational model for cardiomyocytes mechano-electric stimulation to enhance cardiac tissue regeneration. Mathematics8, 1875. 10.3390/math8111875

  • CrossRef
  • Google Scholar

• 35

  Vanegas-AcostaJ. C.Garzón-AlvaradoD. A.ApmZ. (2012). Mathematical model of electrotaxis in osteoblastic cells. Bioelectrochemistry88, 134–143. 10.1016/j.bioelechem.2012.08.002

  • CrossRef
  • Google Scholar

• 36

  WatersC. M.BasslerB. L. (2005). Quorum sensing: Cell-to-cell communication in bacteria. Annu. Rev. Cell Dev. Biol.21, 319–346. 10.1146/annurev.cellbio.21.012704.131001

  • CrossRef
  • Google Scholar

• 37

  WuD.LinF. (2011). A receptor-electromigration-based model for cellular electrotactic sensing and migration. Biochem. Biophysical Res. Commun.411, 695–701. 10.1016/j.bbrc.2011.07.004

  • CrossRef
  • Google Scholar

• 38

  ZajdelT. J.ShimG.WangL.Rossello-MartinezA.CohenD. J. (2020). Scheepdog: Programming electric cues to dynamically herd large-scale cell migration. Cell Syst.10, 506–514.e3. 10.1016/j.cels.2020.05.009

  • CrossRef
  • Google Scholar

• 39

  ZhaoM.PuJ.ForresterJ. V.McCaigC. D. (2002). Membrane lipidsEGF receptors and intracellular signals colocalize and are polarized in epithelial cells moving directionally in a physiological electric field. FASEB J.16, 857–859. 10.1096/fj.01-0811fje

  • CrossRef
  • Google Scholar

• 40

  ZhaoZ.WattC.KarystinouA.RoelofsA. J.McCaigC. D.GibsonI. R.et al (2011). Directed migration of human bone marrow mesenchymal stem cells in a physiological direct current electric field. Eur. Cell. Mat.22, 344–358. 10.22203/eCM.v022a26

  • CrossRef
  • Google Scholar