Numerical Methods to Compute Stresses and Displacements from Cellular Forces: Application to the Contraction of Tissue

We consider a mathematical model for wound contraction, which is based on solving a momentum balance under the assumptions of isotropy, homogeneity, Hooke's Law, infinitesimal strain theory and point forces exerted by cells. However, point forces, described by Dirac Delta distributions lead to a singular solution, which in many cases may cause trouble to finite element methods due to a low degree of regularity. Hence, we consider several alternatives to address point forces, that is, whether to treat the region covered by the cells that exert forces as part of the computational domain or as 'holes' in the computational domain. The formalisms develop into the immersed boundary approach and the 'hole approach', respectively. Consistency between these approaches is verified in a theoretical setting, but also confirmed computationally. However, the 'hole approach' is much more expensive and complicated for its need of mesh adaptation in the case of migrating cells while it increases the numerical accuracy, which makes it hard to adapt to the multi-cell model. Therefore, for multiple cells, we consider the polygon that is used to approximate the boundary of cells that exert contractile forces. It is found that a low degree of polygons, in particular triangular or square shaped cell boundaries, already give acceptable results in engineering precision, so that it is suitable for the situation with a large amount of cells in the computational domain.


Introduction
Wound healing is a complicated, but crucial biological mechanism. In this manuscript, we consider wound healing after skin injury. Since severe (burn) injuries involve a considerable loss of soft tissue, secondary healing takes place. It involves the formation of a blood clot, in case of a cutaneous wound, the regeneration of collagen (extracellular matrix), and revascularisation (which is the re-establishment of a small blood vessel network); see [1] for a biological overview. One of the side effects of secondary healing that follows after a serious skin trauma, is skin contraction. Skin contraction takes place as a result of mechanical, pulling forces that are exerted by the cells (i.e. mainly fibroblasts and myofibroblasts) that are responsible for the regeneration of collagen [2]. Contractions can result in a significant, temporary, or even permanent decrease of area or volume of the damaged tissue. Reductions by 5-10 % of the original wound area have been observed in human skin and in mammalian skin of rodents, even larger reductions have been observed. Such a reduction of skin area or volume leaves residual stresses and strains in the newly repaired skin, as well as in its direct surroundings. This may cause discomfort or even painful sensations to the patient and in extreme cases, contractions may lead to dysfunctionalities of joints. If a contraction is so extreme that the patient develops a disability, then the contraction is referred to as a contracture.
For many of the biological mechanisms that take place during wound healing, mathematical models have been developed. The current manuscript focusses on the formation of a contraction post wounding. Fibroblasts enter the wound site as a result of chemotaxis due to the TGFbeta gradient. Next to the regeneration of collagen, fibroblasts also exert pulling forces to their immediate environment [3]. In some cases, due to being triggered by the high concentration of TGF-beta, fibroblasts differentiate to myofibroblasts, which are known to exert even larger forces than fibroblasts. These larger pulling forces result into the contraction of the tissue around the injury towards the wound centre [4][5][6].
In the literature, several attempts to model the contraction phenomenon can be found [7][8][9][10][11]. The current manuscript focusses on hybrid models for simulating wound contraction in a small scale, where we consider cells as individual entities. We will consider point forces for modelling the balance of momentum, respectively. The modelling framework will entail Dirac Delta functions (distributions), where these pulse-like forces will lead to singularities of the solution in terms of a lower (local) degree of regularity, even such that the solution no longer falls within the finite-element space in which one looks for the solution. Some of the issues have been treated in [12], [13] and [14], regarding well-posedness and finite-element solutions. The treatment of momentum using point forces that we consider in the current paper was developed in [15], [7] and [8].
The quest of several alternative methods is motivated by finding ways to improve accuracy, and by the need of efficiency to simulate the mechanical processes occurring in the skin after a serious (burn) trauma. There are different approaches that treat point forces on the boundary of a cell. One may include the region covered by the cell as part of computational domain. This idea develops into the immersed boundary approach. On the contrary, the 'hole approach', is based on excluding the cell from the computational domain and treat the cell forces as a boundary condition. In this paper, we will focus on the balance of momentum where inertia is neglected and where we assume Hooke's Law to be satisfied. Further, we will use the infinitesimal strain approach. To the best of our knowledge, this paper is the first study that assesses the relation between the 'hole approach' and the immerse boundary approach both analytically and computationally.
The paper is structured as follows. In Section 2, we will discuss the singularity problem occurring in the solution of partial differential equations. Section 3 investigates the 'hole approach' as an alternative to the immersed boundary method, and consistency between these approaches is verified. For a large number of cells in the computational domain, various polygonal approximations of the cell boundary are discussed. In Section 4, we compare the immersed boundary approach to the 'hole approach' and show the results from the polygonal cell approach using various polygonal degrees. Finally some conclusions are presented.

Boundary Value Problems with Point Source
From the definition of the Dirac Delta function, it immediately follows that there is a singularity in the solution to the partial differential equations(PDEs) in some cases. This singularity causes that the solution is irregular and even unbounded if the dimensionality exceeds one. If the PDEs are solved in an infinite domain with Dirac Delta distributions, the solution is known as Green's function. Inspired by this, hereby, we use the Green's function as an intermediate to determine whether there is a singular solution in a given finite domain. In the following contents, we will investigate the solutions in Laplacian equation and elasticity equation respectively.
Proof. Considering Laplacian equation with Dirac Delta function in an infinite region the solution to which is known as the Green's function iŝ Denote v = u −û and then u is extracted as u = v +û. Combining Eq (2.1) and Eq (2.2), a new boundary value problem is derived: The weak form of (BV P 1 ) is for all φ ∈ H 1 (Ω).
Note that the solution of v is classic, which is a sufficient condition that v is in H 1 space. However, the Green's function is not lying in H 1 , since 0∈Ω ∇û 2 dΩ → ∞ regardless of the dimensions d > 1. Since u =û + v, andû / ∈ H 1 (Ω), it immediately follows that u / ∈ H 1 (Ω). To simplify the equation with E = 1 here, the equations above can be combined to Laplacian equation in one dimension: which contains a solution in H 1 (Ω). For dimensions above one, unfortunately, we have found the Green's function in three dimensions in [16]. Therefore, the theorem only states the situation in three dimensions.
Given an open bounded domain 0 ∈ Ω ⊂ R 3 , and the boundary value problem below: where the strain tensor and stress tensor are defined as respectively. Then there does not exist a solution u ∈ H 1 (Ω) such that u can solve (BV P 3 ).
Proof. From [16], the Green's function in three dimensions is where µ and ν is the second Lamé parameter and the Poisson ratio, and i, j present different coordinates. Further, δ ij represents the Kronecker Delta function. The displacement vector of each coordinate can be expressed bŷ (2.10) Thus, similarly as before, letting v = u −û, then the problem becomes in Ω, σ(v) · n + κv = −(σ(n ·û) + κû), on ∂Ω. ∂û i (x) ∂x j 2 is infinite over the domain Ω containing the original point. Here, we will calculate the integral of ∂û x (x) ∂x 2 as an example: Then we rewrite the equation with spherical coordinates as Integrating with respect to r and noting that the inferior of the integral is 0, then where K i (φ, θ), i = 1, 2 is the expression of φ and θ. For other derivative parts, they end up with the same situation in Eq (2.12), that is, for every part of integral 0∈Ω ∇û dΩ, the integral does not exist. Hence, it can be concluded that the Green's function in isotropic open bounded domain is not in H 1 (Ω), which leads to the consequence that the solution to (BV P 3 ), expressed by u = v +û, is not in H 1 (Ω) either.
Remark 2.2. Theorems 1 and 2 can also be proved for the case of homogeneous Dirichlet boundary conditions.

Mathematical Models of Point Forces in Wound Healing 3.1 The Immersed boundary method in R 2
The (myo)fibroblasts exert pulling forces on their immediate surroundings in the extracellular matrix. These forces are directed towards the cell centre and they cause local displacements and deformation of the extracellular matrix. The combination of all these forces cause a net contraction of the tissue around the region, where the fibroblasts are actively exerting forces. The fibroblasts, which are responsible for the regeneration of collagen, enter the wound area after serious trauma due to chemotaxis. Since after restoration of the collagen, the fibroblasts die as a result of apoptosis (programmed cell death), the forces that they exert on their environment disappear. If the deformations are relatively large, then residual stresses remain and permanent displacements remain. Therefore, we consider two types of forces: temporary forces (f t ) and plastic forces (f p ). Here, we will only treat the temporary forces and the way we treat them has been formalized by [15], [7] and [8].
For the temporary force of cell i, the cell boundary Γ i is divided into line segments in the two-dimensional case. We assume that an inward directed force is exerted at the midpoint of every line segment in the normal direction to the line segment. The total force is a linear combination of every force at every segment. Hence, at time t, the total temporary force is expressed by where T N (t) is the number of cells at time t, N i S is the number of line segments of cell i, P (x) is the magnitude of the pulling force exerted at point x per length, n(x) is the unit inward pointing normal vector (towards the cell centre) at position x, x i j (t) is the midpoint on line segment j of cell i at time t and ∆Γ i,j N is the length of line segment j.
Here, x i (t) is a point on the cell boundary of cell i at time t.
The equation for conservation of momentum over the computational domain Ω is given by: In the above equation inertia has been neglected. We treat the computational domain as a continuous linear isotropic elastic domain. Therefore, we use Hooke's Law: where E is the Young's modulus of the domain, ν is Poisson's ratio and is the infinitesimal strain tensor, that is, The above PDE provides a good approximation if the displacements are relatively small. Further, we define the inner product of two second-order n × n tensors (matrices) A and B as follows: where a ij and b ij are the entries of A and B, respectively. On the outer boundary ∂Ω, we use the following Robin boundary condition where κ is a positive constant representing a spring force constant between the domain of computation and its far away surroundings, and u denotes the displacement vector. Note that if κ → ∞, then u → 0 which represents a fixed boundary, and κ → 0 represents a free boundary in the sense that no external force is exerted on the boundary. For the case of only one cell i in the computational domain, we need to solve the following boundary value problem: in Ω, σ · n + κu = 0, on ∂Ω. (3.5) Let V (Ω) be a completion of the Hilbert space H 1 (Ω) with smooth functions [14], then the corresponding weak form of Eq (3.5) on Ω is

The 'Hole Approach' in R 2
Since the force is actually applied on a continuous curve, rather than working on the complete computational domain, we remove the region occupied by the cell. It leaves the computational domain with a hole that is occupied by the cell. Then the force on the cell boundary is modelled by a boundary condition on the boundary of the hole (cell). Therewith, we have boundary conditions on the external boundary, as well as a force boundary condition on the boundary of the cell. The boundary value problem we are working on becomes where n(x) is the unit normal vector pointing out of Ω\Ω C , Ω is the complete computational domain including the cell and extracellular regions, Ω C is the region occupied by the cell, and ∂Ω C is the boundary of the cell. The corresponding weak form for Eq (3.6) is Note that to this problem, it can be proved by combining Lax-Milgram's lemma with Korn's Inequality that a unique solution in H 1 (Ω) exists. In the analysis to come, we assume that the cell stays at the same position and keeps the same shape, hence we have x(t) = x.
Proposition 3.1. Let u H and u I , respectively, be solutions to the 'hole approach' (see Equation (3.6)), and to the immersed boundary approach (see Equation (3.5)). Let ∂Ω C denote the line over which internal forces are exerted, and let ∂Ω be the outer boundary of Ω. Then as ∆Γ −→ 0, Proof. To prove that the above equation holds true, we integrate the PDE of both approaches over the computational domain Ω.
For the immersed boundary approach, we get then after applying Gauss Theorem in the LHS and simplifying the RHS, we obtain By substituting the Robin's boundary condition and sending N i S → ∞, i.e. ∆Γ i,j → 0, the equation becomes Subsequently, we do the same thing for the 'hole approach'. Then, we get − Ω ∇ · σdΩ = 0, and we apply Gauss Theorem: Using the boundary conditions, we get which is exactly the same as Eq (3.7). Hence we proved that Hence, the two different approaches are consistent in the sense of global conservation of momentum and therefore the results from both approaches should be comparable. The only difference between the two approaches is that the 'hole approach' does not consider the stiffness of the cell, since the cell is treated as a hole in the domain. The immersed boundary method contains the internal stiffness of the cell. Therewith, if the cell stiffness is sent to zero, the two formalisms should deliver the same results. Hereby, we are going to prove this transition mathematically and we will see that numerical computations indeed confirm this behaviour.
Before we state and prove a proposition that asserts the transition, we introduce the following energy norm: , where κ is a positive constant. Note that the energy norm is a proper norm according to the definition of norm in [17]. Proposition 3.2. Numerical approximations based on simplicial, continuous finite-element basis functions, to the weak forms of the immersed boundary approach in Equation (3.5) and the 'hole approach' in Equation (3.6), yield the same results upon using the following stiffness for the immersed boundary approach where E is a constant, Ω C is the cell region, Ω\Ω C is the extracellular region and Ω C is surrounded by Ω.
Proof. Due to the symmetry of the tensor (φ), ∀φ, it follows that Hence, rewriting the weak form of the immersed boundary approach taking N i S → ∞, i.e. ∆Γ i,j → 0, (W F I ) becomes

Substituting Eq (3.8) into the above weak form, implies that
Hence, the weak form for the adjusted immersed boundary approach, denoted by (W F I ) is given by: Recalling the weak form of the 'hole approach': We are aware that due to the singularity caused by Dirac Delta distributions in the immersed boundary approach, the solution is no longer in H 1 (Ω). Therefore, following the procedure of discretizing the continuous function space in [12], we approximate the solution by the finite element space V h (Ω) ⊂ H 1 (Ω), such that the solution of (W F I ) can be found in this subset that consists of simplex-based basis functions that are continuous. Subsequently, (W F I ) is given by Applying the same discretizing procedure on the weak form of the 'hole approach', we derive the updated weak form as follows: Note that the above weak forms are identical. Next we demonstrate that the solutions are necessarily the same (hence not determined up to a function or a constant). Since we want to prove the consistency of these two approaches, we rewrite u h in ( Since φ h is a test function, which we can choose freely, such that the provided integrals make sense; Since the energy norm is a proper norm, it can be concluded that v h = 0, in Ω. Hence, we have proved u h I = u h H in Ω. In Proposition 3.2, we have proved the convergence between the finite element solutions to the adjusted immersed boundary approach and the 'hole approach'. Next to it, we are going to prove the convergence between the finite element solution to the adjusted immersed boundary approach and the (exact) solution to the 'hole approach'.

Polygonal Cell Approach
If we consider a domain in which many cells are moving and exerting forces, then the aforementioned two approaches will be very expensive from a computational point of view. Therefore, we will simplify the cell boundary to a low-order polygon, such as to a triangle or square. Furthermore, if the cell size is smaller than the mesh size, then we cannot break the cell boundary into finite segments by the mesh for both approaches. Inspired by finite boundary segments which actually build up a polygon, we can simulate the circular cell by different kinds of polygons.
Eq (3.5) is still used as the basis for the computation of the forces that are exerted by the cells. However, we study the use of just a few boundary segments per cell in such a way that the total force exerted by the cell is the same regardless the order of the polygon.
The cells exert forces on their immediate environment and hence all the points of the computational domain will be displaced. The displacement vector will induce a contraction of the near cell region. This contraction is quantified by the area of the near-cell region. According to [18], for each nodal point, the new position is where X stands for the initial position and x(t) is the position at time t. Defining the gradient matrix of displacement J = ∇ X u, the matrix notation can be worked out as dx = ∂x ∂X dX = (I + ∇ X u)dX = (I + J )dX, (3.10) where ∂x ∂X is the Jacobian matrix. The volume can be calculated by: dx = det(I + J )dX, (3.11) that is, theoretically where Ω 0 is the initial domain. However, to compute the area in Eq (3.12) numerically, we need to apply quadratures like Newton-Côtes quadrature or Gaussian quadrature, which increase the computation expense if we want to track the area at each iteration. Thus, to improve the computational efficiency, another possibility to compute the area of Ω is based on connecting all the nodal points on the boundary to build up a polygon. Then this polygonal area is an approximation of the deformed area since the displacement of each nodal point is available. To calculate the polygon area, one can use shoelace method derived by [19] in 1769. Suppose we have a polygon with n vertices, then the area is calculated by where (x i , y i ), i = 1, . . . , n is the coordinate of vertex i and (x n+1 , y n+1 ) = (x 1 , y 1 ). Note that the vertices should be sorted in counter clockwise or clockwise direction.
To have a better insight of how these different computational approaches affect the cell and the near-cell region, we calculate the reduction of the area with respect to the initial area. If we denote the area after deformation by A Ω and the original area is A 0 Ω , then the ratio is calculated by (3.14) 4 Numerical Results

The Immersed Boundary Approach and The 'Hole Approach'
We use the finite element method to analyse the performance of the immersed boundary approach and 'hole approach'. Since we are interested in the behaviour of the solution in the vicinity of the positions where point forces are exerted, we introduce a subdomain Ω w near the locations where the point sources are exerted. This near-by subdomain, as well as the entire computational domain and the circular line where the forces are exerted are shown in Figure  4.1. The meshes for the two approaches are the same, except for the use of a 'hole' in the hole-approach. The circular curve where the forces are applied models a cell boundary, with its inner region modelling a myofibroblast that exerts forces on its direct environment.
(a) The immersed boundary approach (b) The 'hole approach' The values of the parameters used in this simulation have been listed in Table 4.1. Note that all these parameter values are only for testing the sensitivity of the approaches.
We compare the results from the immersed boundary approach to the results from the 'hole approach'. Figure 4.2 displays the initial cell in blue and the nearby region which is included in the red square, as well as its deformations in black curves. It can be seen that there is a large difference between the results from the two approaches. In particular, the magnitude of the displacement from the 'hole approach' is more than 13 times as large as the displacement from the immersed boundary approach. This discrepancy is caused by the interaction with the region inside the circular cell, which is incorporated in the immersed boundary approach and not in the 'hole approach'. Therefore, we adjust the stiffness of the region inside the circular where γ is a small positive constant. In the following contents about the adjusted immersed boundary approach, we use γ = 10 −5 if there is no further declaration. Then we redo the simulations and plot the results in Figure 4.3. The results of area and total strain energy in the subdomain Ω w have been documented in Table 4.2, and as a result of the use of Eq (3.8), it can be seen that the 'hole approach' and the adjusted immersed boundary approach are consistent since the area reductions are less than a percent. Further, it can be observed that the order of accuracy of the 'hole approach' is slightly better, whereas the adjusted immersed boundary approach is about a factor of four more economical from a computational efficiency point of view.
(a) The immersed boundary approach (b) The 'hole approach'  Due to multiple choices of γ, the value of γ determines the accuracy and convergence of the adjusted immersed boundary approach. In this manuscript, to investigate the effect of γ, it varies from 10 −6 to 10 −3 with steps of a factor of 10. In Table 4.3, besides the area reduction, the convergence rate of the L 2 -norm of the solution and the total strain energy in Ω w are shown. It can be concluded that the value of γ does have a modest impact in the current range, and the influences on various categories are distinct. In other words, for the  Eq(4.1)) and the 'hole approach' (Eq (3.6)) when the same mesh structure used except the hole and the same parameter values applied (Table 4.1). The black line shows the deformed cell and Ω w and the other colour lines represent the original status. area reduction, it is verified that the smaller value γ is, the closer the result is to the one in 'hole approach'. Nevertheless, there is 'bell shape' behaviour appearing for the convergence rate of u L 2 , although the differences are not strikingly large. Further, we observed that, in the perspective of the strain energy in Ω w , the larger γ is, the better the convergence rate.

Polygonal Cell Approach
In the applications that we study, we are interested in multiple cells that are migrating through the computational domain. In typical situations, the cell size is much smaller than the domain size and the cell size could even be smaller than the element size from the discretization. Therefore, it is expensive from a computational point of view to divide the cell boundary into many mesh points and line segments in these applications. Hence, we are interested in the numerical accuracy if each cell is approximated by a simple polygon like a triangle or square instead of a high order polygon. In the presence of multiple small cells, we will study the impact of the polygonal order on the numerical results. The values of the input parameters are given in Table 4.4. In the multi-cell simulations, we locate the cells according to a Point Poisson Process with rate parameter λ, where we choose λ = 15 from [20]. The cell radius has been scaled down to 0.1 of the radius in the previous calculations. The computational domain and the near-cell region are the same as in the earlier simulations. In order to visualize the deformation of the cell and the subdomain Ω w , we set the magnitudes of the forces exerted by the cells to 10. In the simulations, we use the immersed boundary method with low order polygonal approximations of the circular cells. We investigate the performance in terms of the numerical solution with respect to the degree of polygons. An example of a simulation is shown in Figure 4.4, where multiple cells are shown as circles, and the contraction of the region is shown. The cell size is smaller than the mesh size, so we applied the polygonal cell approach here to investigate the area reduction of the region.
The numerical numbers that we investigate are the area reduction due to the pulling forces exerted by the cells and the computation time. In all the calculations where we vary the degree of the polygonal approximation of the cells, we use the same number of cells and the same positions of the centres of the cells. Upon increasing the degree of the polygon, one gradually converges to a circle. In the current computations, we use a maximum number of eight nodes on the cells, that is, we use octagons as the highest polygonal order. The smallest order of polygonal approximation is the triangular shape. We selected the polygons such that the area of each cell is equal in all simulation runs as well as the centres of the cells. Figure 4.5 displays the computation time and relative reduction of area as a function of polygonal degree with multiple cells. Lower order of polygonal approximation admits the advantage that computation time can be reduced due to a lower number of function evaluations from point forces. In the computations, it has turned out that the use of triangles gave a reduction of computation time of roughly fifty percent with respect to the octagonal representation of the cell boundaries according to the histogram in Figure 4.5. The dash line in Figure 4.5 shows that a triangle or square representation of the circles already reproduces the results of the octagonal representation very well, since there is tiny fluctuation. In one word, due to the efficient computation time and good reproduction of the octagonal results in area reduction, we recommend to approximate the cell boundary by a triangle or square if a large number of small cells are used.

Discussion and Conclusions
In this paper, we mainly discussed different approaches to solve linear elasticity problems with point sources forces that are exerted on cell boundaries. In order to simulate wound contraction, it is crucially important to solve the equations for balance of momentum. The body forces are determined by (myo)fibroblasts that exert forces on their immediate extracellular environment. Since we model the forces by the use of point forces which makes the solution not be in the H 1 Sobolev space for dimensions exceeding one, we analysed the relation between the immersed boundary approach and the 'hole approach' and it has been computationally illustrated that the transition from the immerse boundary to the 'hole approach' has a continuous nature with respect to the elasticity in the cellular region. We proved that the finite-element approximations of the two approaches are the same if the stiffness in the cell is neglected. For large numbers of (migrating) cells, it becomes very beneficial to reduce the polygonal order of the representation of the cell boundary. The results indicate that an approximation of a cell boundary by a triangle or square is already sufficiently accurate, and the triangular representation is the least timeconsuming. Furthermore, the computation of the subdomain area by the use of connecting all the boundary vertices to compute a 'polygon' area is the most efficient procedure, combined with applying shoelace method.