The force of cell-cell adhesion in determining the outcome in a nonlocal advection diffusion model of wound healing

1.   Introduction

The movement of cells in spatial micro-environments can be modeled by continuous variables for the continuum density of the cells. The partial differential equations of these models incorporate two components of cell movement: diffusion corresponding to random motion and nonlocal advection corresponding to directed motion. In this work we employ nonlocal advection diffusion models to describe the movement of biological cells in wound healing experiments. In these experiments, cells have a directed motion corresponding to cell-cell adhesion. Cell-cell adhesion is a fundamental feature of multi-cellular organisms, and is linked to maintenance of order in tissue formation, stability, and breakdown.

There are several signaling pathways that control normal cell processes, such as cell proliferation, cell division, cellular adhesion, and apoptosis. One of the most important signaling pathways is transforming growth factor beta $TGF\beta$, which regulates the proliferation, differentiation, adhesion, migration, and apoptosis of many cell types, including endothelial cells, hematopoietic cells, lymphocytes, and ECM production.

In wound healing experiments, cells lie in a viscous thin layer in a laboratory well, and are subjected to experimental drug therapies. Cells in the layer are moving and proliferating. The velocity of a cell is proportional to the net adhesive force acting on it, due to bonds formed with nearby cells. This adhesive flux results, because cells sense their surroundings over a sensing radius via protrusions such as filopodia, which interact with nearby cells to facilitate intercellular communication. The sensing radius is considerable larger than the cell radius. The division of cells depends on a nutrient furnished to the cell population. In the experiments, the layer of cells is scored laterally to form a narrow depression, which then restores completely to the original level layer, or incompletely to a series of ridges.

In our references we have cited works that treat laboratory wound healing experiments ^{[1,2,3,4,5,6,7,8,9,10]}. In Figure 1, images from wound healing experiments in ^[4] are given. In these experiments, human umbilical vein endothelial cells were plated on gelatin-coated plastic dishes and wounded with a pipette tip, and then imaged by microscope over several hours.

In ^[11], numerical simulations were given to investigate the dependence of wound closure on the diffusion parameter $\alpha$. In this work, our objective is to examine the relationship in (2.1), between the sensing radius $\rho$ and the adhesive force parameter $\lambda$ in $g(n) = n \, (\lambda - n)$. We will demonstrate numerically, that this relationship depends on $\rho$ and $\lambda$ in such a way that for $\lambda$ sufficiently large, there is a critical value $\rho_c(\lambda)$ of $\rho$ such that $n(x, t)$ converges to $\bar{n}(x)$ for $\rho$ less than $\rho_c(\lambda)$. For $\rho$ greater than $\rho_c(\lambda)$, $n(x, t)$ does not converge to $\bar{n}(x)$, and incomplete closure occurs, with a series of ridges established above and below $\bar{n}(x)$. We will also examine the spacing $\rho_s(\lambda)$ of the spatial ridges in the case of incomplete wound closure, as a function of $\lambda$.

2.   Nonlocal advection diffusion model applied to wound healing experiments

In our references we have cited works that treat laboratory wound healing experiments with mathematical models ^{[2,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,
51,52,53,54,55,56,57,58,59,60,61,62,63,64,65]}. We note in particular the nonlocal advection diffusion models of wound healing experiments by Stephen Gourley ^[11,25,56]. In ^[11], the partial differential equation of the nonlocal advection diffusion model for wound healing has the following form:

where $x$ is spatial position, $-\infty < x < \infty$, and $t$ is time, $t \geq 0$.

The model in ^[17] for cell-cell adhesion has a direct representation of cell-cell contact using a nonlocal spatial term in the model equation. This formulation was originally due to Armstrong, Painter, and Sherratt ^[13], and subsequently applied to somite formation in early vertebrate development (Armstrong, Painter, and Sherratt ^[14]), to tumour growth (Gerisch and Chaplain ^[39]), and to cancer invasion (Sherratt, Gourley, Armstrong, and Painter ^[56]). These authors considered basic properties of the models, such as positivity, boundedness, and spatial pattern formation, which are important for biological applications.

In (2.1), $n(x, t)$ is the 1-dimensional continuum density of cells with respect to the lateral position $x$ across the wound at time $t$. The initial condition $n(x, 0) = n_0(x)$ corresponds to the wound depth at spatial position $x$ at time $0$. The domain of $n(x, t)$ is taken as $(-\infty, \infty)$, since the boundary of the laboratory well is not significant compared to the sensing radius. The interpretation of the model outcomes relates to early or intermediate dynamics over time, which are not influenced by cells near the boundary of the laboratory well.

In (2.1), $\alpha$ is the diffusion constant, $g$ describes the advection force on a cell at position $x$ due to other cells within a sensing radius $\rho$ centered at $x$, $h$ describes variation of the advection force over the sensing radius centered at $x$, and $f$ describes the loss or gain of cells at position $x$. The sensing radius quantifies the total forces attributable to bonds formed by cells at nearby locations due to cell protrusions such as filopodia. The magnitude of these local forces depends on the number of attachments made by cells within the experimentally specified sensing radius.

In (2.1), $\int_{-\rho}^{\rho} \, g(n(x+\hat{x}, t) \, h(\hat{x}) \, d \hat{x}$ represents the total adhesive force, attributable to cell-cell bonding, on a cell at position $x$ from cells at position $x + \hat{x}$, $\hat{x} \in (x - \rho, x + \rho)$. The function $g(n)$ is approximately linear for small $n$, decreases for $n$ above some threshold, and is effectively $0$ above some critical value. In our numerical simulations of (2.1), we will assume that $g(n) = n (\lambda - n)$ with $\lambda > 1$. We will call $\lambda$ the cell-cell adhesion constant.

The function $h(\hat{x})$ in (2.1) represents the magnitude and direction of the adhesive force over the sensing radius of the cell as a function of $\hat{x} \in (-\rho, \rho)$. This force diminishes with increasing separation of cells. The adhesive force is directed towards cell centers. The force on a cell at $x$ due to a cell at $x+\xi$ will tend to pull the cell at $x$ to the right if $\xi > 0$, and to the left if $\xi < 0$. The direction of the adhesive force, toward the cell at position $x$, requires $h(\hat{x})$ to be an odd function, with $\hat{x} \geq 0$ for $\hat{x} \geq 0$. In our numerical simulations, we will take $h(\hat{x}) = 0.1 \arctan(0.2 \, \hat{x})/\arctan(0.2)$. The function $f(n)$ in (2.1) is of logistic type, so that at lower cell densities, there is increased cell division, and at higher cell densities there is cell loss due to overcrowding. In our simulations, we will take $f(n) = n (1 - n)$, corresponding to a normalization of the density at the equilibrium value $\bar{n}(x) \equiv 1$.

In ^[11], the existence, regularity, positivity, boundedness, and asymptotic behavior of the solutions of (2.1) was investigated with the above assumptions. In ^[25], the model was extended to higher dimensional spatial domains for the cell densities. In ^[30], the model was extended to include proliferating cell age structure corresponding to the cell cycle. In ^[30], operator semigroup methods were used to establish suffucient conditions for the local asymptotic stability of the wound free equilibrium. In ^[11], the following sufficiency result was proven for the global asymptotic stability of the equilibrium $\bar{n}(x) \equiv 1$ in the model (2.1):

Theorem. Let $h$ be an odd function such that $h(\xi)\geq 0$ when $\xi\geq 0$, let $f(n) = n(1-n)$, and let $g(n) = n(\lambda-n)$ with $\lambda > 1$. Suppose that $\delta < n_0(x) \leq \lambda$ for all $x \in R$, where $\delta > 0$. Suppose

Then, $n(x, t)$ converges to $1$ as $t \rightarrow \infty$, exponentially in $L^2(R)$.

In general, the result is sufficient, but not necessary for complete would healing.

3.   Numerical simulations of the wound healing experiments

In our numerical simulations we will view cells in the laboratory well as a culture grown to confluence in nutrient medium, and then incubated in serum-free medium. Before the wound is imposed, the cell density is constant at the normalized equilibrium value $\bar{n}(x) \equiv 1$. The layer of cells is scraped with a plastic pipette tip, and then replenished with serum-free medium. The healing of the wound can be viewed as the lateral closure of the wound across the scoring. Dependent on model parameters, the proliferating cells close the wound back to the equilibrium density $\bar{n}(x)$, or form an incomplete closure consisting of a series of ridges above and below the equilibrium.

In our simulations, the diffusion parameter $\alpha = 1.0$, the initial condition is taken as $n_0(x) = 1.0 - 0.75 \, e^{- (0.1 x)^6}$, corresponding to the lateral depression of the cells forming the wound, centered at $x = 0$ (Figure 2), and the adhesion force magnitude function is $h(\hat{x}) = 0.1\, \arctan(0.2 \, \hat{x}) / \arctan(2.0)$. Comparable results can be obtained for similar values of $\alpha$, $n_0(x)$, and $h(\hat{x})$.

With the parameters and initial condition above, we will illustrate wound closure dependence on the sensing radius $\rho$ and the adhesion force parameter $\lambda$ in the force of adhesion function $g(n) = n \, (\lambda - n)$. We will determine the critical value $\rho_c(\lambda)$ distinguishing complete wound closure from incomplete wound closure. We will examine the character of the spatial density function $n(x, t)$ in the case of incomplete wound closure with respect to the spacing of the oscillating spatial ridges of $n(x, t)$. We will provide three examples of $\lambda$ to illustrate the critical sensing radius $\rho_c(\lambda)$ as a function of $\lambda$: $\lambda = 5.0$, $\lambda = 6.0$, and $\lambda = 7.0$.

In the first simulation, $\lambda = 5.0$ and $\rho = 16.9$. The graph of $n(x, t)$ is given for time from $t = 0$ to $t = 300$ in Figure 3. In Figure 4, the $n(x, t)$ is graphed at $t = 2, t = 10, t = 70, t = 300$. As time advances, $n(x, t)$ converges to the equilibrium $\bar{n}(x) \equiv 1.0$, and the wound closes completely. At every time level, the graph is concave up at $x = 0.0$.

In the second simulation, $\lambda = 5.0$ and $\rho = 17.0$. The graph of $n(x, t)$ is given for $t = 0$ to $t = 200$ in Figure 5. In Figure 6, $n(x, t)$ is graphed at $t = 2, \, t = 10, t = 70, t = 200$. As time advances, $n(x, t)$ does not converge to the equilibrium $\bar{n}(x) \equiv 1.0$, and the wound does not close completely. In the early dynamics, the wound stabilizes to a series of ridges with spacing $\approx 38.5$. At every time level, the graph of $n(x, t)$ is concave up at $x = 0.0$. The critical value $\rho_c(\lambda)$ for $\lambda = 5.0$ is $\rho_c(5.0) \approx 17.0$. The three sufficiency conditions above for $\lim_{t \rightarrow \infty} n(x, t) = 1$, require $\rho < 1.72$ with $\delta = 0.25 = n_0(0)$.

In the third simulation, $\lambda = 6.0$ and $\rho = 10.9$. The graph of $n(x, t)$ is given for $t = 0$ to $t = 300$ in Figure 7. In Figure 8, $n(x, t)$ is graphed at $t = 2.0, \, t = 10.0, t = 70.0, t = 300.0$. As time advances, $n(x, t)$ converges to the equilibrium $\bar{n}(x) \equiv 1.0$, and the wound closes completely. At every positive time level, the graph is concave down at $x = 0.0$.

In the fourth simulation, $\lambda = 6.0$ and $\rho = 11.0$. The graph of $n(x, t)$ is given for $t = 0$ to $t = 300$ in Figure 9. In Figure 10, $n(x, t)$ is graphed at $t = 2.0, \, t = 10.0, t = 70.0, t = 300.0$. As time advances, $n(x, t)$ does not converge to the equilibrium $\bar{n}(x) \equiv 1.0$, and the wound does not close completely. In the early dynamics, the wound stabilizes to a series of ridges with spacing $\approx 27.8$. At every time level, the graph is concave down at $x = 0.0$. The critical value $\rho_c(\lambda)$ for $\lambda = 6.0$ is $\rho_c(6.0) \approx 11.0$. The three sufficiency conditions above for $\lim_{t \rightarrow \infty} n(x, t) = 1$, require $\rho < 1.18$ with $\delta = 0.25 = n_0(0)$.

In the fifth simulation, $\lambda = 7.0$ and $\rho = 7.9$. The graph of $n(x, t)$ is given for time from $t = 0$ to $t = 100$ in Figure 11. In Figure 12, $n(x, t)$ is graphed at $t = 2.0, \, t = 10.0, t = 70.0, t = 100.0$. As time advances, $n(x, t)$ converges to the equilibrium $\bar{n}(x) \equiv 1.0$, and the wound closes completely. At every time level, the graph is concave down at $x = 0.0$.The critical value $\rho_c(\lambda)$ for $\lambda = 7.0$ is $\rho_c(7.0) \approx 8.0$. The three sufficiency conditions above for $\lim_{t \rightarrow \infty} n(x, t) = 1$, require $\rho < 1.9$ with $\delta = 0.25 = n_0(0)$.

In the sixth simulation, $\lambda = 7.0$ and $\rho = 8.0$. The graph of $n(x, t)$ is given for time from $t = 0$ to $t = 100$ in Figure 13. In Figure 14, $n(x, t)$ is graphed at $t = 2.0, \, t = 10.0, t = 70.0, t = 100.0$. As time advances, $n(x, t)$ does not converge to the equilibrium $\bar{n}(x) \equiv 1.0$, and the wound does not close completely. In the early dynamics, the wound stabilizes to a series of ridges with spacing $\approx 20.0$. At every time level, the graph is concave down at $x = 0.0$. The critical value $\rho_c(\lambda)$ for $\lambda = 7.0$ is $\rho_c(7.0) \approx 8.0$. The three sufficiency conditions above for $\lim_{t \rightarrow \infty} n(x, t) = 1$, require $\rho < 0.86$ with $\delta = 0.25 = n_0(0)$.

In the Table below, we provide the values of the critical sensing radius $\rho_c(\lambda)$ as a function of $\lambda$. We also provide, as a function of $\lambda$, the spacing of the maxima (and minima) $\rho_s(\lambda)$ of the density functions as they oscillate, at the critical sensing radius values $\rho_c(\lambda)$. Both $\rho_c(\lambda)$ and $\rho_s(\lambda)$ are decreasing functions of the force of adhesion parameter $\lambda$. In Figure 15 we graph $\rho_c(\lambda)$ and $\rho_s(\lambda)$.

4.   Conclusions

A model of wound healing experiments has been investigated, based on the model (2.1) in ^[11]. A theoretical analysis was given in ^[11] for the existence, regularity, and boundedness of solutions of (2.1). Additionally, a sufficient condition on the parameters was established, for the convergence of solutions $n(x, t)$ to the complete wound free equilibrium $\bar{n}(x) \equiv 1.0$, as time advances.

In this work, the connection between the force of adhesion $g(n) = n (\lambda - n)$ and the sensing radius of cells $\rho$, is investigated for this model. Numerical simulations were given to connect the sensing radius $\rho$ to the force of adhesion parameter $\lambda$, with the model parameters $\alpha = 1.0$ and initial condition $n_0(x) = 1.0 - 0.75 \, e^{-(0.1 x)^6}$. These numerical simulations demonstrate that for each $\lambda\, \in [4.5, 10.0]$, there is a critical value $\rho_c(\lambda)$ of the sensing radius, such that complete wound closure occurs for $\rho < \rho_c(\lambda)$ and incomplete wound closure occurs for $\rho > \rho_c(\lambda)$. This critical value $\rho_c(\lambda)$ is graphed in Figure 15 and decreases as $\lambda$ increases. For $\lambda \leq 4.0$, the numerical simulations demonstrate that the wound closes for every value of the sensing radius $\rho$.

The result in ^[11] providing sufficient conditions on the sensing radius $\rho$ to yield complete wound healing gives values considerably lower than the critical sensing radius $\rho_c(\lambda)$ for the numerical simulations for $\lambda = 5.0, 6.0, 7.0$. A complete theoretical establishment of necessary and sufficient conditions for the asymptotic stability of the equilibrium $n(x, t) \equiv 1, \, t \geq 0$, corresponding to complete wound closure in (2.1), remains an open problem.

The simulations also demonstrate the character of the ridges formed in incomplete closure, as time advances, with respect to their oscillation spacing. This oscillation spacing value $\rho_s(\lambda)$ is also graphed in Figure 15 and decreases as $\lambda$ increases. These numerical simulations demonstrate the value of the cell–cell adhesion model (2.1) for understanding the output of wound healing experiments, which are of major importance in medical research of therapies for human diseases. The numerical simulations were carried out using MATHEMATICA and are available upon request.

Conflict of interest

The author declares there is no conflicts of interest.