Mathematical model for promotion of wound closure with ATP release

Kenta Odagiri; Hiroshi Fujisaki; Hiroya Takada; Rei Ogawa

doi:10.2142/biophysico.bppb-v20.0023

Abstract

To computationally investigate the recent experimental finding such that extracellular ATP release caused by exogeneous mechanical forces promote wound closure, we introduce a mathematical model, the Cellular Potts Model (CPM), which is a popular discretized model on a lattice, where the movement of a “cell” is determined by a Monte Carlo procedure. In the experiment, it was observed that there is mechanosensitive ATP release from the leading cells facing the wound gap and the subsequent extracellular Ca²⁺ influx. To model these phenomena, the Reaction-Diffusion equations for extracellular ATP and intracellular Ca²⁺ concentrations are adopted and combined with CPM, where we also add a polarity term because the cell migration is enhanced in the case of ATP release. From the numerical simulations using this hybrid model, we discuss effects of the collective cell migration due to the ATP release and the Ca²⁺ influx caused by the mechanical forces and the consequent promotion of wound closure.

Significance

Mechanical forces to cells can cause various biological functions. One interesting example of them is the promotion of wound closure by mechanical stimulation. Extracellular ATP release caused by external stimulation triggers Ca²⁺ influx into the cells, resulting in accelerated wound closure. We here proposed the mathematical model for wound closure with mechanosensitive ATP release and subsequent Ca²⁺ response. Numerical simulations using our model showed that collective cell migration caused by ATP release promotes fast wound closure.

Introduction

Recent studies in the field of mechanobiology [1–3] have revealed that cells recognize and respond to various external stimuli or external forces, and it has been found that various biological functions in vivo are triggered by this mechanobiological effects [4–9]. As such, mechanobiological issues are important and challenging not only for basic biology but also for the optimal treatment of plastic surgery which is called “mechanotherapy” [10–14].

In relation to wound closure and mechanobiology, Takada, Furuya and Sokabe found that the mechanical forces on cells are important for wound closure [15]. In their experiment, a linear wound was made in cultured vascular endothelial cells, and the wound closure process was examined by applying stretch stimulation. It was found that the wound closure rate increased with the stretch stimulation compared to no external stimulation. They also measured real-time imaging of ATP and Ca²⁺ concentrations, where the stretch stimulation releases large amounts of ATP from the cells closest to the wound gap. Here ATP functions as an extracellular messenger and increases intracellular Ca²⁺ levels. When ATP travels like a transient wave from the proximal wound to the posterior cell, it is followed by a transient wave of Ca²⁺. It was shown that the response of the cells to Ca²⁺ was indispensable for the wound closure because the Ca²⁺ wave disappeared and cell migration did not occur when the extracellular Ca²⁺ was removed.

Motivated by this experimental study, we here propose a mathematical model of wound closure process that explicitly considers ATP release by mechanical stimulation to the cells, the accompanying change in Ca²⁺ concentration, and collective cell migration by response to Ca²⁺. We construct a hybrid model of the Cellular Potts Model (CPM) [16–25] and the Reaction-Diffusion (RD) equations [26,27], where CPM represents the dynamics of cells, such as cell migration and cell growth, and RD represents the changes in ATP and Ca²⁺ concentration. The cell migration in response to Ca²⁺ is expressed as self-propulsion driven by cell polarity, which is further added as an energy term in the CPM. Using this hybrid model, we try to numerically simulate the experiment on wound closure by Takada and coworkers. We here focus on the relationship between the ATP release and the rate of wound closure rather than explicitly considering the mechanical force applied to the cells, because the ATP release by mechanical stimulation occurred after a very brief mechanical stimulation.

This paper is organized as follows. In Section “Materials and Methods”, we first introduce the hybrid mathematical model for the wound closure process. We next show the numerical results of our model in Section “Results and Discussion”. To compare with the experimental results such that the wound closure is promoted by ATP release due to the mechanical stimulation, we here examine two different scenarios with and without ATP release. We also discuss why the wound closure rates are different in these two cases. Finally, we conclude this paper in Section “Conclusions”.

Materials and Methods

We present a mathematical model for cellular dynamics in wound closure that incorporates extracellular ATP release and subsequent extracellular Ca²⁺ influx due to mechanical stimulation. Our approach is a hybrid model that couples cellular dynamics with spatiotemporal dynamics of chemical components both inside and outside of the cells. Specifically, we propose a model that combines the Cellular Potts Model with Reaction-Diffusion Equations. The Cellular Potts Model is used to represent cellular population dynamics, while the Reaction-Diffusion Equations capture the spatiotemporal dynamics of ATP and Ca²⁺. We provide further details below.

It is worth noting that some previous studies have incorporated the mechanical interaction with the substrate by coupling the Cellular Potts Model with the Finite Element Method (FEM) [22–24]. However, our study focuses on the effects of ATP release on wound closure and thus does not explicitly consider mechanical stimulation.

Cellular Dynamics

Collective cell migration is one of the key processes in the wound closure process, and various mathematical models have been proposed [28–33]. We here employ the Cellular Potts Model (CPM) for the computational modeling of the cellular collective movement. CPM is defined on a lattice, representing a cell shape as a group of sites on the lattice. Each site has a cell index σ, where the same cell shares the same cell index σ. To mimic various cellular behavior such as cell deformation and migration, CPM repeats the “exchange” of adjacent lattice sites with a Metropolis algorithm using a kind of “energy” as described below.

The move is described as an attempt to copy the cell index σx→ at a randomly selected site x→ into another randomly selected neighboring site x'→. Whether the move is accepted or not is determined using the change in the total energy ΔE due to the exchange. When ΔE≤0, the move is always accepted, however, when ΔE>0, the move is accepted with the following probability: px→→x'→=e-β∆E, where β denotes the magnitude of cell fluctuation [17], in accord with the conventional Metropolis algorithm.

Energy of Cell Population

The energy of cell population used in the CPM calculations is modelled as follows:

ECPM=∑x→, x'→Jτσx→, τσ x'→1-δσx→, σ x'→ +λV∑kVk-Vt,k2+λS∑kSk-St,k2-Emig ,

(1)

where the first term represents the interfacial energy due to the cell adhesion, x→ and x'→ denote two neighboring lattice sites, σx→ denotes a cell index at a lattice site x→, and τ denotes a cell type (Note that the term “cell index” refers to an individual cell, while “cell type” refers to all cells that share the same characteristics). J_τ,τ’ denotes the adhesion energy between cell type τ and τ’. The second term represents cell volume constraint energy, where V_k and V_t,k denote a calculated volume and a target volume of cell k, respectively, and λ_V is the strength of the volume constraint. The third term represents cell surface constraint energy, where S_k and S_t,k denote a calculated surface area and a target surface area of cell k, respectively, and λ_S is the strength of the surface constraint (In two-dimensional system, cell surface area means the length of cell’s perimeter). The second and third terms have the effect of maintaining the proper volume (V_t,k) and surface area (S_t,k) of cell k, respectively. The last term represents a cell migration energy due to cellular self-propulsion as described below.

Cellular Self-Propulsion

The cell migration energy E_mig with the self-propulsion is given by

Emig=∑kPktpk→pk→∙vk→ ,

(2)

where P_k(t) is the strength of the cell autonomous motility of cell k at time t, pk→ represents the cell polarity vector of cell k. vk→ represents the displacement vector of the center of cell k before and after a single site exchange, and is calculated as vk→=gk,af→-gk,be→, where gk,af→ and gk,be→ are the center of cell k after and before single site exchange, respectively. This formulation indicates that the cell polarity drives cellular self-propulsion. We assume that Ca²⁺ ions enhance the cellular self-propulsion, and thus set the cell autonomous motility P_k(t) to depend on the intracellular Ca²⁺ concentration in each cell as Pkt=θCk¯t-CthP in this study. θ(a) is the Heaviside function (1 when a>0 and 0 otherwise), Ck¯t is the average concentration of Ca²⁺ inside cell k at time t, C_th is the threshold value of Ca²⁺ concentration, and P is the magnitude of cell autonomous motility. If P_k(t) is a function that switches motility depending on the intracellular Ca²⁺ concentration, our simulation results would not qualitatively change.

Szabó and Czirók proposed that the cell polarity vector pk→ decays spontaneously but reinforced by cell displacements [34,35]. The variation of pk→ (∆pk→) is calculated based on the following equation at the end of each Monte Carlo step (MCS).

∆pk→=-γpk→+∆gk→ .

(3)

The first term indicates spontaneous decay of pk→ with decay rate γ. The second term indicates reinforcement of pk→ resulting from cell displacements, where gk→ and ∆gk→ are the center of cell k and the change in gk→ during a single MCS.

Dynamics of Chemical Components

The spatiotemporal changes in extracellular ATP concentration Ax→,t and intracellular Ca²⁺ concentration Cx→,t are modelled by the following Reaction-Diffusion equations (RD).

∂Ax→,t∂t=δσx→, σEdgeθtθtb-tk0-δσx→, σcellskAAx→,t+DA∇2Ax→,t ,

(4)

∂Cx→,t∂t=δσx→, σcellsθAx→,t-AthkCC0-Cx→,t .

(5)

The first and second terms on the right-hand side of Equation (4) denote ATP release and ATP consumption. ATP molecules release from only the leading cells σ_Edge that faces the wound edge at the beginning of the cell migration 0≤t≤t_b) and are consumed by cells σ_cells with rate k_A as indicated by the second term. t_b is the time when the ATP release ends and k₀ is the release rate of ATP. The last term in Eq. (4) denotes the extracellular diffusion of ATP with the diffusion coefficient D_A. Equation (5) represents the intracellular Ca²⁺ response in the cells σ_cells induced by released extracellular ATP molecules, where A_th is the threshold value of the ATP concentration for the extracellular Ca²⁺ influx due to the Ca²⁺ response, k_C is the influx rate of Ca²⁺, and C₀ is the steady-state concentration of extracellular Ca²⁺. We numerically solve the above reaction-diffusion equations using the Crank-Nicolson scheme on a lattice that corresponds to the CPM lattice, employing 10 diffusion steps per MCS with time step Δt=0.1 and spatial step size Δx=1.0. The Neumann boundary condition is applied.

Cell Growth and Contact Inhibition

We further introduce the cell growth considering the contact inhibition. Contact inhibition is the mechanism which normal cells stop proliferation or growth when they contact each other.

Cell growth is modelled by increasing the target volume V_t,k and the target surface area S_t,k [36,37] in the CPM. If cell k remains in a state where it is not completely surrounded by other cells for a certain period of time, T_g, the value of target volume V_t,k increase by one. Here, T_g corresponds to the reciprocal of the growth rate of cell volume. Once V_t,k increases by one, an additional time interval of T_g is required until the next increase. V_t,k can increase up to its maximum value V_max. On the other hand, if cell k is completely surrounded by the other cells, the value of V_t,k does not change, indicating that the cell does not get large due to the effect of the contact inhibition. To ensure that the value of V_t,k is increased at the appropriate time during CPM calculations, a timer for cell growth, t_g,k, is set. If cell k is not completely surrounded by other cells and V_t,k is not at its maximum value V_max, the value of timer t_g,k is increased by one each MCS. When the value of timer t_g,k reaches T_g, the value of V_t,k is increased by one and the value of timer t_g,k is reset to zero to measure the time until the next increase in V_t,k. If a cell is completely surrounded by other cells or if V_t,k is at its maximum value, cell growth stops, and the value of timer t_g,k remains unchanged.

Furthermore, assuming that cells tend to adopt an elongated shape during migration, we control the shape of a cell by adjusting the value of the target surface area, S_t,k, based on the value of V_t,k. To achieve an elongated shape, we assume a cell aspect ratio of 4:1 and compute the value of S_t,k from the value of V_t,k. Assuming the shape of a cell to be a rectangle, with the length of the longer side denoted as 4r and the length of the shorter side denoted as r, the two-dimensional volume (area) V can be calculated as V=4r², and the two-dimensional surface area (perimeter) S can be calculated as S=10r. By considering these V and S values as V_t,k and S_t,k, respectively, the following relationship can be derived: St,k=5Vt,k. When the value of V_t,k is updated, the value of S_t,k is also updated based on this relationship.

It should be noted that this model does not take into account cell division, as the time scale of the experimental results [15] is shorter than the time scale of cell division (cells do not divide during the simulation time).

Coupling the Cellular Potts Model with the Reaction-Diffusion Equations

We here describe the method used to couple CPM representing cell dynamics with RD representing chemical dynamics for computation. Since the time scale of chemical dynamics is much shorter than that of cell dynamics, we execute the calculation of RD ten times during a single Monte Carlo step (MCS) of the CPM. The summary of the computational procedure for each MCS is presented in the following pseudocode algorithm.

1. Perform the following site-exchange calculations for the number of times equal to the number of lattice sites (system size) in the CPM.

a. Randomly select of two sites for site-exchange.

b. Calculate the change in the total energy ΔE for site-exchange using the equations (1) and (2).

c. Perform the Metropolis algorithm using the calculated ΔE.

2. Numerically calculate the RD equations (4) and (5) ten times.

3. Perform the following cell growth calculations for each cell.

FOR k=1 to Number of Cells DO

a. Calculate the equation (3) for updating the cell polarity vector pk→

b. IF (timer for cell growth t_g,k≥period T_g) THEN

IF (V_t,k<V_max) AND (cell k is not completely surrounded by other cells) THEN

V_t,k=V_t,k+1 and then St,k=5Vt,k (update target volume V_t,k and target surface area S_t,k)

t_g,k=0 (reset timer t_g,k to zero to measure the time until the next increase in V_t,k)

END IF

ELSE (timer for cell growth t_g,k<period T_g) THEN

IF (V_t,k<V_max) AND (cell k is not completely surrounded by other cells) THEN

t_g,k=t_g,k+1 (increment the timer for cell growth t_g,k)

END IF

END DO

Results and Discussion

We here show some numerical results using the hybrid model (CPM with polarity+RD equations for chemicals+contact inhibition for cells) for the wound closure process. To mimic the experimental situation [15], we set up the initial distribution of cells as shown in Figure 1 (a). The gap region (grey) between the upper and lower cells represents the linear wound. Green cells in Figure 1 (a) indicates the leading cells that face the wound edge. Only these leading cells release ATP molecules, as in the experiment [15], at the beginning of the simulation (cell migration) until t_b=1 (we estimate 1 MCS as 15 sec). The initial volume of each cell is 25 (5×5) units (20×20 μm²), where we consider a square lattice with 1×1 unit (4×4 μm²), and the initial number of the cells is 288. The total system size is 120×120 units (480×480 μm²), and the wound region is the half of it (480×240 μm²). This wound gap length corresponds to the experimental study [15].

Figure 1

(a) Initial distribution of cells. F, M and B represent cells facing the wound gap, cells located in the middle of the cell cluster and cells located at the back of the cell cluster, respectively. (b) and (c) Time evolution of wound closure process in the ATP-free model and the ATP model, respectively.

Table 1 shows the parameters of the numerical simulation (units are arbitrary). J_CC and J_CM denote the cell-cell adhesion energy and the cell-medium (cell free area) adhesion energy, respectively. We set the energy related parameters, such as J_CC, J_CM, λ_V, λ_S, P, and γ, to stabilize cell shape and allow deformation, growth, and migration to occur within the designated simulation time frame. Additionally, the ratio between the magnitude of cell fluctuation β and the change in the total energy ΔE is also crucial for cell deformation, migration, and stability of cell shape. Hence, we set β=0.2 in our simulations to enable a certain level of cell deformation and motility, while also ensuring the stability of cell shape. We set the values of the parameters related RD equations (4) and (5), taking into consideration that ATP and Ca²⁺ rapidly spread within a few seconds in the experiments [15]. We assumed that it takes approximately 1.5 hours (360 MCS) at the shortest to reach the maximum size (V_max) in this simulation since we have no experimental data on the time required for cell growth. In this simulation, we set T_g=10 as the time required for a cell to increase its size by one. As a result, the time required for a cell to grow from its initial size of 25 to its maximum size of 60 is 350 MCS (350=10×(60–25)), which is in consistent with our previous assumption of 360 MCS.

Table 1 Parameters used in the simulations (arbitrary units)

Symbol	Meaning of parameter	Value	Equation to be used
β	magnitude of cell fluctuation	0.2	Metropolis calculation
J_CC	cell-cell adhesion energy	9.0	Eq. (1)
J_CM	cell-medium adhesion energy	3.0	Eq. (1)
λ_V	strength of the volume constraint	2.0	Eq. (1)
λ_S	strength of the surface constraint	0.5	Eq. (1)
P	magnitude of cell autonomous motility	2.0	Eq. (2)
C_th	threshold value of [Ca²⁺] for cell autonomous motility	0.1	Eq. (2)
γ	decay rate of cell polarity vector pk→	0.025	Eq. (3)
t_b	time when the ATP release ends	1	Eq. (4)
k₀	ATP release rate from cells	10.0	Eq. (4)
k_A	ATP consumption rate of cells	0.01	Eq. (4)
D_A	diffusion coefficient of ATP	0.1	Eq. (4)
A_th	threshold value of ATP concentration for Ca²⁺ influx	0.1	Eq. (5)
k_C	Ca²⁺ influx rate	1.0	Eq. (5)
C₀	steady-state concentration of extracellular Ca²⁺	1.0	Eq. (5)
T_g	time until the next increase in the value of target volume	10	Cell growth calculation
V_max	maximum value of target cell volume	60	Cell growth calculation

To clarify the effect of the collective cell migration due to the ATP release and the subsequent Ca²⁺ influx resulting from the mechanical forces, we here numerically examine two different scenarios with and without the ATP release: The former is referred to as ATP model and the latter as ATP-free model. In the ATP-free model, the ATP release is inhibited by setting t_b=0. As a result, Ca²⁺ response does not occur, leading to E_mig=0 in the equation (2). Compared to our previous study of the wound closure [38], we add to this model the effect of contact inhibition and investigate in detail the effects of the ATP release. Therefore, this study is more suitable for qualitative comparison with the previous experimental study [15].

Wound Closure by the Collective Cell Migration

Figure 1 (b) and (c) show the time evolution of the wound closure process in the ATP-free model and ATP model, respectively. These figures clearly show that the wound gap closes more quickly in the ATP model than in the ATP-free model. In the ATP model, the cells facing the wound gap rapidly migrate forward, and this quick migration causes rapid closure of the wound gap. On the other hand, in the ATP-free model, the leading cells migrate forward very slowly, and thus the wound gap closes very slowly.

To compare the rate of wound closure in these models, we measure the change in the surface coverage by cells over time (Figure 2). Although the rate of wound closure in the two cases is almost the same at the early phase (t=0~200), the rate in the ATP-free model suddenly slows down after around t=400. On the other hand, in the ATP model, the wound closure rate does not slow down significantly until cells fill most of the wound gap, and the wound gap almost fill (Surface coverage reaches 0.90) at around t=1000 (4.2 hour). Therefore, the ATP release clearly shows the promotion of the wound closure with the parameters we chose in Table 1.

Figure 2

Time evolution of the surface coverage by cells in the ATP-free model and the ATP model.

ATP Diffusion and Ca²⁺ Response

We next show the ATP diffusion and the subsequent Ca²⁺ response at the early stage in the ATP model. In the previous experimental study [15], the diffusion of released ATP and the subsequent Ca²⁺ response have been observed within a few seconds after mechanical stimulation. Figure 3 (a), (b) and (c) show the time evolution of the distribution of cells, ATP and Ca²⁺ at the early stage, respectively. The released ATP from the leading cells diffuses quickly and the resulting Ca²⁺ influx occurs although the width of the wound gap is almost the same (t=10~50). And after a while, the diffused ATP decays (t=100,250). These results qualitatively correspond to the experimental results.

Figure 3

Time evolution of the distribution of cells (a), ATP (b) and Ca²⁺ (c) at the early stage. Red (white) indicates that the concentration of the chemical components is high (low), and gray indicates that the concentration of them is zero (almost zero). ATP diffusion and the subsequent Ca²⁺ influx occur quickly although the wound gap is barely closed.

Differences in the Wound Closure Rate

To investigate the difference in the wound closure rate, we here show the degree of cell migration at various location within the cell cluster and the distribution of the cell size.

We first calculate the average cell position (perpendicular to the wound gap) at various location within the cell cluster, that is, cells in the row facing the wound gap (cluster F in Figure 1 (a)), cells located in the middle row within the cell cluster (cluster M in Figure 1 (a)) and cells located in the back row within the cell cluster (cluster B in Figure 1 (a)), respectively. Figure 4 (a), (b) and (c) show the time evolution of the average cell position in each cluster in two models. Time evolution of the average position in cluster F is quite similar to that of the surface coverage by cells (Figure 2). It indicates that the cell migration of the leading cells directly reflects the wound closure rate. On the other hand, the average cell positions in clusters M and B in the two models show completely different behaviors. In the ATP model, cells in clusters M and B also migrate forward, whereas they remain almost the same positions in the ATP-free model. It shows that quick cell migration for cluster F due to the ATP release causes the subsequent cell migration for clusters M and B in the ATP model.

Figure 4

(a), (b) and (c) correspond to time evolution of the average cell position in clusters F, M and B, respectively.

In Figure 1, it may appear that each cell moves separately to close the wound gap. However, in reality, the cells move collectively in an aligned orientation. To quantitatively verify the collective cell migration due to the ATP release, we introduce the following local polar order parameter, S_Ψ(t) for a cell cluster Ψ [4],

SΨt=1NΨ∑i∈Ψμi→μi→ ,

(6)

where N_Ψ is the number of cells in cluster Ψ and μi→ is the velocity vector of i-th cell in cluster Ψ. Figure 5 shows the time evolution of S_Ψ averaged over a moving average for 10 intervals. In cluster F, the value of S_Ψ in the ATP model is greater than that in the ATP-free model until t=750 (Figure 5 (a)). It indicates that cells in cluster F in the ATP model move in a more aligned orientation. Therefore, the ATP release promotes the collective cell migration for cluster F.

Figure 5

(a), (b) and (b) Time evolution of the local polar order parameter in clusters F, M and B, respectively.

We next visualize the distribution of the cell size in each cell cluster. Figure 6 (a) and (b) is basically the same as Figure 1 (b) and (c) but the color code is different: Differences in cell color indicate the differences in cell size. Blue color means the size of a cell is large (V≃V_max), and white color means the size of a cell is small (V≤25). Figure 6 (a) clearly indicates that only the cells in cluster F grow and the size of the cells in cluster M and B (far from the wound) hardly changes in the ATP-free model. To investigate it quantitatively, we calculate the average cell size at various locations in clusters F, M and B. Although cells in cluster F in two models have almost the same size over time (Figure 7 (a)), cells at M and B in the ATP model grow faster than those in the ATP-free model (Figure 7 (b) and (c)). This is the reason why the increase in the surface coverage slows down in the ATP-free model.

Figure 6

(a) and (b) correspond to Figure 1 (a) and (b), respectively, and blue color indicates large cells whereas white color does small cells.

Figure 7

(a), (b) and (c) Time evolution of the average cell size in clusters F, M and B, respectively.

We next investigate why cells inside the cluster in the ATP-free model can hardly grow. In our numerical model, a cell surrounded by other cells cannot get large due to the contact inhibition. Therefore, vacant ratio of a cell, which denote how much space is available around the cell, is a useful for determining whether the cell can grow or not. Figure 8 (a), (b) and (c) show the time evolution of the average vacant ratio in clusters F, M and B, respectively. In the ATP-free model, cells inside the cluster (clusters M and B) are almost completely surrounded by other cells, resulting in an average vacant ratio that is almost zero or very small. As a result, the cells cannot grow large due to contact inhibition. It should be noted that the decrease in the average vacant ratio in cluster F in the ATP model after t=700 is due to the rapid migration and growth of the leading cells, which results in the loss of vacant space.

Figure 8

(a), (b) and (c) Time evolution of the average vacant ratio in clusters F, M and B, respectively.

From the above considerations, we summarize the reason for the promotion of wound closure in the ATP model as follows. In the ATP model, the value of E_mig in the equation (2) becomes significantly large because of the high intracellular Ca²⁺ concentration, making it easier for the front cells to move and for space to spread rapidly, thus enabling the rear cells to migrate and grow smoothly. In contrast, in the ATP-free model, E_mig becomes 0 due to the lack of Ca²⁺ in the cells, causing the front cells to move less and thus preventing the rear cells from migrating or growing. Therefore, the collective cell migration induced by the ATP release and the Ca²⁺ influx promotes the cell growth inside the cell cluster and the consequent fast wound closure.

Conclusions

We have introduced a hybrid model which consists of the CPM (cell dynamics) and the RD equations (dynamics for chemical components) for the wound closure process, and further added the cell polarity and contact inhibition. We here focused on the wound closure rate for the ATP-free model (without ATP release from the cells facing the wound) and the ATP model. The difference in the wound closure rate is caused by the rapid cell migration due to the ATP release and the subsequent Ca²⁺ influx induced by the mechanical forces, which we model as the polarity of the cell movement. The rapid cell migration makes space for the cells inside the cell cluster, and consequently causes rapid growth for these cells. In contrast, in the ATP-free model, the leading cells migrate slowly, and thus the cells inside the cluster do not have enough space to get large and cannot grow due to the contact inhibition.

Our hybrid model can also be applied to the vascular network formation, where anisotropic shape of a cell is important to generate such a fibrous form. By incorporating the effect of cell shape and taking into account the effect of the ATP release induced by external forces, our hybrid model will be further extended to deal with the vascular network formation under mechanobiological conditions. And also, it would be essential to add the effect of mechanical interaction with cells and the substrate to our hybrid model for more realistic simulations.

In this study, many model parameters were set based on assumptions due to the difficulty in correlating them with experimental results. A future task is to estimate the model parameters using data assimilation methods based on large amounts of cell imaging data, similar to those used in the field of meteorology.

Conflict of Interest

The authors declare no conflict of interest.

Author Contributions

Conceptualization, K.O, H.F.; methodology, K.O., H.F., H.T., R.O.; validation, K.O., H.F., H.T., R.O.; investigation, K.O., H.F.; data creation, K.O.; writing---original draft preparation, K.O.; writing---review and editing, K.O., H.F., H.T., R.O. All authors have read and agreed to the published version of the manuscript.

Data Availability

The evidence data generated and/or analyzed during the current study are available from the corresponding author on reasonable request.

A preliminary version of this work, DOI: https://doi.org/10.48550/arXiv.2209.01354, was deposited in the arXiv on September 3, 2022.

Acknowledgements

We are grateful to Kazue Kudo, Yoshitaro Tanaka, and Shigeyuki Komura for inspiring and useful discussions. This research was partially supported by the Japan Society for the Promotion of Science (KAKEN 22K11941 to H.F., K.O., 19K12207 to K.O., H.F., H.T.) and AMED-CREST, AMED (JP20gm0810012) to K.O., H.F., H.T., R.O.

References

[1] Jacobs, C. R., Huang, H., Kwon, R. Y. Introduction to cell mechanics and mechanobiology (Garland Science, New York, 2012).
[2] De, S., Hwang, W., Kuhl, E. Multiscale modeling in biomechanics and mechanobiology (Springer London, 2015).
[3] Rawlinson, S. C. Mechanobiology: Exploitation for Medical Benefit (John Wiley & Sons, Hoboken, 2017).
[4] Hayakawa, K., Tatsumi, H., Sokabe, M. Actin stress fibers transmit and focus force to activate mechanosensitive channels. J. Cell Sci. 121, 496–503 (2008). https://doi.org/10.1242/jcs.022053
[5] Ingber, D. E. Tensegrity-based mechanosensing from macro to micro. Prog. Biophys. Mol. Biol. 97, 163–179 (2008). https://doi.org/10.1016/j.pbiomolbio.2008.02.005
[6] Jansen, K. A., Donato, D. M., Balcioglu, H. E., Schmidt, T., Danen, E. H., Koenderink, G. H. A guide to mechanobiology: Where biology and physics meet. Biochim. Biophys. Acta Mol. Cell Res. 1853, 3043–3052 (2015). https://doi.org/10.1016/j.bbamcr.2015.05.007
[7] Ladoux, B., Mège, R. M. Mechanobiology of collective cell behaviours. Nat. Rev. Mol. Cell Biol. 18, 743–757 (2017). https://doi.org/10.1038/nrm.2017.98
[8] Huang, G. Mechanobiology in wound healing. Biophys. J. 121, 173–174 (2022). https://doi.org/10.1016/j.bpj.2021.12.016
[9] He, S., Li, X., Ji, B. Mechanical force drives the polarization and orientation of cells. Acta Mech. Sin. 35, 275–288 (2019). https://doi.org/10.1007/s10409-019-00864-z
[10] Thompson, W. R., Scott, A., Loghmani, M. T., Ward, S. R., Warden, S. J. Understanding mechanobiology: Physical therapists as a force in mechanotherapy and musculoskeletal regenerative rehabilitation. Phys. Ther. 96, 560–569 (2016). https://doi.org/10.2522/ptj.20150224
[11] Akaishi, S., Akimoto, M., Ogawa, R., Hyakusoku, H. The relationship between keloid growth pattern and stretching tension: visual analysis using the finite element method. Ann. Plast. Surg. 60, 445–451 (2008). https://doi.org/10.1097/SAP.0b013e3181238dd7
[12] Huang, C., Holfeld, J., Schaden, W., Orgill, D., Ogawa, R. Mechanotherapy: Revisiting physical therapy and recruiting mechanobiology for a new era in medicine. Trends. Mol. Med. 19, 555–564 (2013). https://doi.org/10.1016/j.molmed.2013.05.005
[13] Huang, C., Du, Y., Ogawa, R. Mechanobiology and mechanotherapy for cutaneous wound‐healing. in Mechanobiology: Exploitation for Medical Benefit. pp. 239–253 (John Wiley & Sons, Hoboken, 2017). https://doi.org/10.1002/9781118966174.ch15
[14] Fu, S., Panayi, A., Fan, J., Mayer, H. F., Daya, M., Khouri, R. K., et al. Mechanotransduction in wound healing: From the cellular and molecular level to the clinic. Adv. Skin Wound Care 34, 67–74 (2021). https://doi.org/10.1097/01.ASW.0000725220.92976.a7
[15] Takada, H., Furuya, K., Sokabe, M. Mechanosensitive ATP release from hemichannels and Ca2+ influx through TRPC6 accelerate wound closure in keratinocytes. J. Cell Sci. 127, 4159–4171 (2014). https://doi.org/10.1242/jcs.147314
[16] Graner, F., Glazier, J. A. Simulation of biological cell sorting using a two-dimensional extended Potts model. Phys. Rev. Lett. 69, 2013–2016 (1992). https://doi.org/10.1103/PhysRevLett.69.2013
[17] Glazier, J. A., Graner, F. Simulation of the differential adhesion driven rearrangement of biological cells. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 47, 2128–2154 (1993). https://doi.org/10.1103/physreve.47.2128
[18] Merks, R. M., Glazier, J. A. A cell-centered approach to developmental biology. Physica A 352, 113–130 (2005). https://doi.org/10.1016/j.physa.2004.12.028
[19] Hirashima, T., Hosokawa, Y., Iino, T., Nagayama, M. On fundamental cellular processes for emergence of collective epithelial movement. Biol. Open 2, 660–666 (2013). https://doi.org/10.1242/bio.20134523
[20] Guisoni, N., Mazzitello, K. I., Diambra, L. Modeling active cell movement with the Potts model. Front. Phys. 6, 61 (2018). https://doi.org/10.3389/fphy.2018.00061
[21] Czirók, A., Varga, K., Méhes, E., Szabó, A. Collective cell streams in epithelial monolayers depend on cell adhesion. New J. Phys. 15, 075006 (2013). https://doi.org/10.1088/1367-2630/15/7/075006
[22] van Oers, R. F., Rens, E. G., LaValley, D. J., Reinhart-King, C. A., Merks, R. M. Mechanical cell-matrix feedback explains pairwise and collective endothelial cell behavior in vitro. PLoS Comput. Biol. 10, e1003774 (2014). https://doi.org/10.1371/journal.pcbi.1003774
[23] Rens, E. G., Merks, R. M. Cell contractility facilitates alignment of cells and tissues to static uniaxial stretch. Biophys. J. 112, 755–766 (2017). https://doi.org/10.1016/j.bpj.2016.12.012
[24] Rens, E. G., Merks, R. M. Cell shape and durotaxis explained from cell-extracellular matrix forces and focal adhesion dynamics. iScience 23, 101488 (2020). https://doi.org/10.1016/j.isci.2020.101488
[25] Rens, E. G., Edelstein-Keshet, L. From energy to cellular forces in the Cellular Potts Model: An algorithmic approach. PLoS Comput. Biol. 15, e1007459 (2019). https://doi.org/10.1371/journal.pcbi.1007459
[26] Kapral, R., Showalter, K. Chemical waves and patterns (Kluwer, Dordrecht, 1995).
[27] Kuramoto, Y. Chemical Oscillations, Waves, and Turbulence (Springer-Verlag, Berlin, 1984 (Dover Edition, 2003)).
[28] Palachanis, D., Szabó, A., Merks, R. M. Particle-based simulation of ellipse-shaped particle aggregation as a model for vascular network formation. Comput. Part. Mech. 2, 371–379 (2015). https://doi.org/10.1007/s40571-015-0064-5
[29] Suzuki, T. Mathematical methods for cancer evolution (Springer, Singapore, 2017).
[30] Barton, D. L., Henkes, S., Weijer, C. J., Sknepnek, R. Active vertex model for cell-resolution description of epithelial tissue mechanics. PLoS Comput. Biol. 13, e1005569 (2017). https://doi.org/10.1371/journal.pcbi.1005569
[31] Ishihara, S., Marcq, P., Sugimura, K. From cells to tissue: A continuum model of epithelial mechanics. Phys. Rev. E 96, 022418 (2017). https://doi.org/10.1103/PhysRevE.96.022418
[32] Méhes, E., Vicsek, T. Collective motion of cells: from experiments to models. Integr. Biol. (Camb) 6, 831–854 (2014). https://doi.org/10.1039/c4ib00115j
[33] Alert, R., Trepat, X. Physical models of collective cell migration. Annu. Rev. Condens. Matter Phys. 11, 77–101 (2020). https://doi.org/10.1146/annurev-conmatphys-031218-013516
[34] Szabó, A., Ünnep, R., Méhes, E., Twal, W. O., Argraves, W. S., Cao, Y., et al. Collective cell motion in endothelial monolayers. Phys. Biol. 7, 046007 (2010). https://doi.org/10.1088/1478-3975/7/4/046007
[35] Giverso, C., Scianna, M., Preziosi, L., Buono, N. L., Funaro, A. Individual cell-based model for in-vitro mesothelial invasion of ovarian cancer. Math. Model. Nat. Phenom. 5, 203–223 (2010). https://doi.org/10.1051/mmnp/20105109
[36] Stott, E. L., Britton, N. F., Glazier, J. A., Zajac, M. Stochastic simulation of benign avascular tumour growth using the Potts model. Math. Comput. Model. 30, 183–198 (1999). https://doi.org/10.1016/S0895-7177(99)00156-9
[37] Szabó, A., Merks, R. M. Cellular potts modeling of tumor growth, tumor invasion, and tumor evolution. Front. Oncol. 3, 87 (2013). https://doi.org/10.3389/fonc.2013.00087
[38] Odagiri, K., Fujisaki, H. Mathematical model for wound healing caused by exogeneous mechanical forces. AIP Conf. Proc. 2343, 020017 (2021). https://doi.org/10.1063/5.0048360

責任著者(Corresponding author)

J-STAGEへの登録はこちら（無料）