Point Forces in Elasticity Equation and Their Alternatives in Multi Dimensions

We consider several mathematical issues regarding models that simulate forces exerted by cells. Since the size of cells is much smaller than the size of the domain of computation, one often considers point forces, modelled by Dirac Delta distributions on boundary segments of cells. In the current paper, we treat forces that are directed normal to the cell boundary and that are directed toward the cell centre. Since it can be shown that there exists no smooth solution, at least not in $H^1$ for solutions to the governing momentum balance equation, we analyse the convergence and quality of the approximation. Furthermore, the expected finite element problems that we get necessitate scrutinizing alternative model formulations, such as the use of smoothed Dirac distributions, or the so-called smoothed particle approach as well as the so-called 'hole' approach where cellular forces are modelled through the use of (natural) boundary conditions. In this paper, we investigate and attempt to quantify the conditions for consistency between the various approaches. This has resulted in error analyses in the $H^1$-norm of the numerical solution based on Galerkin principles that entail Lagrangian basis functions. The paper also addresses well-posedness in terms of existence and uniqueness. The current analysis has been performed for the linear steady-state (hence neglecting inertia and damping) momentum equations under the assumption of Hooke's law.


Introduction
Wound healing is a complicated process of a sequence of cellular events contributing to resurfacing, reconstitution and restoration of the tensile strength of injured skin.Significant damage of dermal tissue often leads to skin contraction.If the contraction of the skin near a joint is large then it may result into a decrease of functionality.If the patient's daily life is impacted as result of the contraction, then one speaks of a contracture.
In order to improve the patient's quality of life, one aims at reducing the contractile behavior of the skin.To reduce the severity of the contraction, one needs to know the physiological dynamics and time evolution of the underlying biological mechanisms.According to [5,7,10], the contraction starts developing during the proliferative phase of wound healing.This proliferative phase sets in after the inflammatory phase, in which the immune system is clearing up the debris that resulted from the damage.The proliferative phase usually starts from the second day post-wounding, and commonly lasts two to four weeks.Besides the closure of the epidermis (that is the top layer of skin), the proliferative phase is characterized by ingress of fibroblasts from the surrounding undamaged tissue and differentiation to myofibroblasts, and by the regeneration of collagen by the (myo)fibroblasts.Despite the relatively quick closure of the epidermis, often the restoration of the underlying dermis is still in progress.After closure of the epidermis, the damaged region in the dermis is referred to as a scar instead of a wound.Next to the regeneration of collagen, the (myo)fibroblasts exert contractile forces on their direct surroundings, which will result into contraction of the scar tissue.In human skin, typically volume reductions of 5 -10% are commonly observed [6].
The current manuscript contains an extension of the work in [8], which treats a model for the contractile forces exerted by the (myo)fibroblasts.The forces are distinguished into two categories: (1) temporary forces that are exerted as long as the (myo)fibroblasts are actively pulling; and (2) permanent or plastic forces, which are imaginary forces that are introduced to describe the localized plastic deformations of the tissue.This formalism was firstly developed by Vermolen and Gefen [14], and later extended by Boon et al. [3].The formalism is based on the point forces, which are mathematically incorporated by means of linear combinations of Dirac Delta distributions.The irregular nature of Dirac Delta distributions make the solution to the elliptic boundary value problem from the balance of momentum have a singular solution in the sense that for dimensionality higher than one, no formal solutions in the finite-element space H 1 exist.Although in classical finite-element strategies, one uses for instance piecewise linear Lagrangian elements, of which the basis functions are in H 1 , and therewith one attempts to approximate the solution (which is not in H 1 ) as well as possibly by a function in H 1 .Bertoluzza et al. [2] demonstrated the convergence of finite-element solutions by means of piecewise linear Lagrangian elements in multiple dimensions.In our earlier studies [11,12], we proved the convergence of solutions obtained by regularization of Dirac Delta distributions, the so-called smoothed particle approach and the so-called 'hole' approach to the solution obtained by Dirac Delta distributions in the one-and two-dimensional cases.In the one-dimensional case, for the sake of completeness, we start with the presentation of force equilibrium with point forces, the equations are given by To simplify the equation, we use E = 1 here, the equations above can be combined to the one-dimensional Laplace equation: We assume that there is a biological cell with size h and centre position c in the computational domain such that 0 < c − h/2 < c < c + h/2 < L. Then the force is given by f Combined with homogeneous Dirichlet boundary conditions: the Galerkin form is given by where (x) + = max{0, x}.Note that in one dimension, the solution is piecewise linear and hence in H 1 (Ω), however not in H 2 (Ω).Since most conventional errors are expressed in the L 2 -norm of the second derivative of the solution, one may not apriorily expect very accurate finite element solutions.
In the current manuscript we extend the results to general dimensionality.The boundary value problem is stated in Section 2. The 'hole' approach and the smoothed particle approach are developed in Section 3. Furthermore, we prove consistency between all the alternatives and the immersed boundary approach in multi dimensions.Section 5 displays some conclusions and discussions.

Elasticity Equation with Point Sources in Multi Dimensions
Let Ω be a bounded domain in R n , then we consider the following balance of momentum where inertial effects have been neglected: Here σ denotes the stress tensor and f represents a body force that is exerted within Ω.We consider a linear, homogeneous, isotropic and continuous material; hence, Hooke's Law is used here for the relation between the stress and strain tensors: where E is the stiffness of the computational domain, ν is Poisson's ratio and is the infinitesimal Eulerian strain tensor: Within the domain of computation, Ω, we consider the presence of a biological cell, which occupies the portion Ω C that is completely embedded within Ω (hence Ω C is a strict subset of Ω).The boundary of the cell Γ C is divided into surface elements.On the centre of each surface element, a point force by means of Dirac Delta distributions, is exerted in the direction of the normal vector that is directed inward into the cell.This results into (see [14]): where N S is the number of surface elements of the cell, P (x, t) is the magnitude of the pulling force exerted at point x and time t per unit of measure (being area in R 3 or length in R 2 ), n(x) is the unit inward pointing normal vector (towards the cell centre) at position x, x j (t) is the midpoint on surface element j of the cell at time t and ∆S(x j ) is the measure of the surface element j.In the general model where we use this principle, we consider transient effects due to migration and possible deformation of the cells.However, since we predominantly focus on the mathematical issues in the current manuscript, we will not consider any time-dependencies and hence t will be removed from the expressions in the remainder of the paper.
In the n-dimensional case, we are solving the boundary value problems described in Eq (2.1), (2.2) and (2.3).The body force is given in Eq (2.4).Therefore, the immersed boundary value problem that we are going to consider is given by Next to this boundary value problem, we consider the continuous immersed boundary counterpart, given by where we take N s → ∞.Thus, the body force reads as Due to the irregular nature of the Dirac Delta distributions, the solutions do not exist in H 1 .We attempt to approximate the solution by the functions in H 1 via the Galerkin form of (BV P ) and (BV P ∞ ).In this manuscript, piecewise linear Lagrangian basis functions are selected.Further, the convergence of the finite-element solutions using linear Lagrangian elements in general dimensionality has been proved in [2].
To construct the Galerkin form, we introduce the bilinear form a(., .) where the last step is motivated by symmetry of the stress tensor σ.Since the solution u is not in H 1 (Ω), we consider a subspace of H 1 (Ω), which is defined as V h (Ω) = Span{φ 1 , φ 2 , . . ., φ N } [13].Here, φ i for i = {1, 2, . . ., N } is the linear Lagrangian basis function, which is piecewise smooth and continuous over Ω.Hence, these basis functions are in H 1 .Subsequently, the Galerkin form is We further consider the solution to the continuous immerse boundary problem, with the following Galerkin form: Before we proceed to claim the existence and the uniqueness of the Galerkin solution in (GF ), we state Korn's Inequality in multiple dimensions: Lemma 2.1.(Korn's Second Inequality [4]) Let Ω ⊂ R n be an open, bounded and connected domain.Then there exists a positive constant K, such that for any vector-valued function u ∈ H 1 0 (Ω), We note that Korn's Second Inequality can be generalised to cases in which the boundary condition u = 0 is imposed only on a non-zero measure part of the boundary.Using Korn's Second inequality gives the following lemma: Lemma 2.2.Let Ω ⊂ R n be an open, bounded and connected domain.Then there exists a positive constant K, such that for any vector-valued function u ∈ H 1 0 (Ω), Proof.The lemma directly follows from the definition of the stress tensor, let u ∈ H 1 0 (Ω): The last step follows from Lemma 2.1.Hence, redefining K := E 1+ν K concludes the proof the lemma.
Herewith, coerciveness of the linear form a(., .)has been demonstrated, which is needed for the proof of existence and uniqueness of the Galerkin finite-element solution.
Theorem 2.1.Let {φ i } be piecewise Lagrangian basis field functions and let F be a vector in R n with unit length, further let P ∈ C(Ω), and let j=1 P (x j )n(x j )φ h (x j )∆S(x j ) for all φ h ∈ V h , and u h = N S j=1 P (x j )u G h (x; x j ; n(x j ))∆S(x j ); Proof.
• It is immediately clear that a(., .) is a bilinear form.We have V h ⊂ H 1 0 (Ω), and a(., .) is bounded in H 1 0 (Ω) (see for instance [1]).Furthermore, Lemma 2.2 says that a(., .) is coercive in H 1 0 (Ω).Regarding the right-hand side, we have |φ h | ≤ M 1 for some M 1 > 0 since φ h is a Lagrangian function, and hence the magnitude of the right-hand side can be bounded from above by Note that ||F || = 1.Hence the right-hand side is bounded, since we are looking for a solution in a finite dimensional space V h , the system where the coefficients of the symmetrix matrix A are defined by a ij = a(φ i , φ j ), and where a limited number of entries of b are non-zero and given by F • φ h (x ), which is finite.Since b is finite, and A is invertible, existence and uniqueness of u h follow (one could apply Lax-Milgram's theorem on the space R n in this context) from the algebraic system.
• Existence and uniqueness follow analogously, only boundedness of the right-hand side, which is a linear functional in φ h ∈ V h (Ω) has to be checked: Note that n has unit length.The summation gives the polygonal length or polyhedral area of the cell boundary.Hence the right-hand side is bounded, then by Lax-Milgram's Lemma, existence and uniqueness follow.Further by substitution, it follows that that The last step uses the first part of the theorem, and finally the assertion is proved similarly to the first assertion.
• We proceed similarly, by boundedness of the right-hand side: where |Γ C | is the measure of the boundary surface of the biological cell.It again shows that the right-hand side is a bounded linear functional in V h (Ω).We proceed by substitution: Note that, formally, it was not necessary to prove boundedness, since coerciveness implies uniqueness and the existence was proved by construction and by combining the result for the existence of u G h .Note that for the 'continuous' weak formulation, there is no solution in H 1 , hence the above claim demonstrates the existence and uniqueness of a Galerkin-based approximation in a subset of H 1 to a function that is not in H 1 .The situation is somewhat comparable to approximating √ 2 / ∈ Q arbitrarily accurately by a sequence of successive approximations in Q.Further in twoand three-dimensional case, the convergence between the solution to (GF ) and (GF ∞ ) can be proved.Similar work has been done in [9] regarding Stokes problem with the Delta distribution term.
Theorem 2.2.Let Γ C be a polygon or polyhedron embedded in Ω ⊂ R n and let P (x) be sufficiently smooth.Further, let x j be the midpoint of surface element ∆S(x j ).Denote u ∆S h as the Galerkin solution to (GF ) and the u ∞ h as the Galerkin solution to (GF ∞ ), respectively.In two dimensions, for any x / ∈ Γ C , there exists a positive constant C, such that for each component of u ∞ h we have , where h max is the maximal diameter among all the triangular elements over Γ C .
Proof.Away from Γ C , the function u G h is smooth, and since P (x) is smooth as well, the integrand, given by P (x)u G h is smooth as well.For ease of notation, we set f (x) = P (x)u G h (x; x ; n).We start with the 2D-case.Given the i − th boundary element ∆S i on Γ C with the endpoints x i and x i+1 and we denote its midpoint by x i+1/2 , where i ∈ {1, 2, • • • , N S }.We consider , and subsequently We calculate the contribution over ∆S i to the integral, where Taylor's Theorem and the Mean Value Theorem for integration are used to warrant the existence of a ŝ ∈ (−1, 1), such that where H(x(s)) is the Hessian matrix of f (x(s)).Therefore, we obtain that Since f (x) ∈ C 2 (Ω), it follows that there exists a K > 0, such that Therefore, considering the summation of the boundary elements over ∂Ω C , where ∆S max = max i∈{1,...,N S } ||x i+1 − x i || is the maximal length of the line segment over Γ C , and |Γ C | is the perimeter of the polygon Γ C .It can be concluded that there exists a positive constant K, such that In three dimensions, the surface element is a triangle.We map the triangle in (x, y, z)space to the reference triangle in (s, t)-space with points (0, 0), (0, 1) and (1, 0).Suppose there is a surface element e j with nodal points x 1 , x 2 and x 3 , then the centre point of e j is x c = (x 1 + x 2 + x 3 )/3.The map from the reference triangle e 0 to the physical triangle e j is given by For any function f (x) ∈ C 2 (Ω), the integral over the original triangle is given by where J is the Jacobian matrix, given by and | det(J T J )| is twice the area of the original triangle e j , i.e.
We conduct the same process as for the two dimensional case, we obtain, where x( 1 3 , 1 3 ) = x c coincides with the midpoint of element e j , and where Taylor's Theorem for multi-variate functions is used: Due to f (x) ∈ C 2 (Ω), then for the Hessian matrix of f (x), there exists K > 0, such that It yields where h 2 max is the largest diameter in the original triangle e j .Considering all the surface elements over Γ C , we compute where h 2 max is the maximal diameter among all the surface element (i.e.triangle) and |Γ C | is the sum of the measure (area in R 3 ) of all the surface elements over Γ C .Therefore, in three dimensions, we can conclude that there exists a positive constant K, such that for the unique Galerkin solution to both (GF ) and (GF ∞ ), The above proof and theorem can easily be extended to higher dimensionalities.

Alternative Approaches for Elasticity Equation with
Point Sources in Multi Dimensions

The 'Hole' Approach
A different approach is based on considering cellular forces on the cell boundary by means of a boundary condition.In this alternative approach, one 'removes' the cell region from the domain of computation.Herewith, one creates a 'hole' in the domain.We consider the balance of momentum over Ω \ Ω C .This gives the following boundary value problem: where σ is defined in Eq (2.2) with stiffness E. Let D ⊂ Ω, then we introduce the following notation: Note that the stiffness can be a constant or a function of space over the domain D.
The corresponding weak form is stated below: Since φ ∈ H 1 (Ω\Ω C ), it follows from the Trace Theorem [4], and by noting that φ| ∂Ω = 0, that there is a , which implies that the right-hand side in the weak form is bounded.Subsequently one combines Korn's Inequality with Lax-Milgram's Lemma to conclude that a unique solution in H 1 exists.We compare the immersed boundary method with the 'hole' approach by taking β 0, then we adjust the immersed boundary method such that Regarding the adjusted immersed boundary approach where the stiffness is given by Eq (3.1), we have the following Galerkin form where V h (Ω) is defined in Theorem 2.1 in Section 2.
For the 'hole' approach, we have the following Galerkin form We will prove that the adjusted immersed boundary method is a perturbation of the 'hole' approach: Proposition 3.1.Let u H h and u β h , respectively, satisfy Galerkin forms (GF H ) and (GF β ), then there is a Proof.First we note that, as in the spirit of Theorem 2.1, we consider Galerkin solutions in a subset of H 1 whereas the solution to the 'continuous' weak formulation is not in H 1 .Formally (GF H ) and (GF β ) hold for test functions φ h from different sets, namely V h (Ω) and V h (Ω \ Ω C ).If we choose V h (Ω C ) to correspond to Lagrangian basis functions associated to internal nodes in Ω C , then these basis functions vanish at Γ C .Furthermore, within the set of Lagrangian basis functions that are associated with Ω \ Ω C , there are Lagrangian basis functions associated with Γ C , which have a compact, hence limited, support over Ω C and in Ω\Ω C , then let v = u β h −u H h , then subtraction of problems (GF H ) and (GF β ) gives The left-hand side is a bounded and coercive form on which we can apply Korn's Inequality.Furthermore, boundedness of the right-hand side in V h (Ω \ Ω C ) follows by application of the Cauchy-Schwartz Inequality, hence there is an . Herewith, we arrive at Note that the a Ω\Ω C (v, φ h ) contains v and φ h in Ω \ Ω C , whereas the right-hand side of the inequality contains norms over Ω C .Using Korn's Inequality, and upon setting φ h = v in Ω\Ω C , we arrive at For the case of a spring-force boundary condition on ∂Ω one can derive a compatibility condition.To this extent, we consider the following boundary value problems, for the 'hole' problem: and for the immersed boundary problem: Next we give a proposition regarding compatibility for the 'hole' approach and the immersed boundary method for the case of a spring boundary condition: Let u H and u I , respectively, be solutions to the 'hole' approach, see (BV P H ) and to the immersed boundary approach, see (BV P I ).Let Γ C denote the boundary of the cell, over which internal forces are exerted, and let ∂Ω be the outer boundary of Ω.Then 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 letting N S → ∞, i.e. ∆S(x j ) → 0, the equation becomes Subsequently, we do the same thing for the 'hole' approach.Then, we get and we apply Gauss Theorem: Using the boundary conditions, we get which is exactly the same as Eq (3. 2).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 Smoothed Particle Approach
The Gaussian distribution is used here as an approximation for the Dirac Delta distribution.Hereby, we show that in the n-dimensional case, the Gaussian distribution is a proper approximation for the Dirac Delta distribution.

Lemma 3.1. For an open domain
Firstly, we integrate over the infinite domain: Again let s i = , and furthermore , i = {1, 2, . . ., n}.We denote x 1 = (x 1,1 , x 2,1 , . . ., x n,1 ), x 2 = (x 1,2 , x 2,2 , . . ., x n,2 ) and x = (x 1 , x 2 . . ., x n ).By Taylor Expansion, f (x) can be rewritten as where H(x ) is Hessian matrix of f (x).For any non-negative integer d, x n,1 x n,1 x n,1 s n,1 π exp(−y) for y > 0 and the fact that exp(y) < 1 y α as y → ∞, we see that the second term approximates zero faster than the first term.Hence, we conclude that there is a C > 0 such that As a remark we add that setting f (x) = 1, immediately shows that there is a Using the result above, we start with analysing different approaches with only one relatively big cell in the computational domain.According to the model described in Eq (2.4), the forces released on the boundary of the cell are the superposition of point forces on the midpoint of each surface element.For example, if we use a square shape to approximate the biological cell, then the forces are depicted in Figure 3.1.Therefore, in n dimensional case (n > 1), if the biological cell is a n-dimensional hypercube, then the forces can be rewritten as where e i is the standard basis vector with 1 in the i-th coordinate and 0 s elsewhere, and ∆x is the length of cell boundary in each coordinate.For the smoothed force approach, we set δ(x) ≈ δ ε (x).The force is given by (3.4) Following the same process in two dimensions [12] and thanks to the continuity of Gaussian distribution, as ∆x → 0, the force converges to f S = P (∆x) n ∇δ ε (x − x ). (3.5) and u ε h be the Galerkin solution to for all φ h ∈ V h (Ω). (3.7) Then there is an Using bilinearity of a(., .)gives upon setting w = u h − u ε h the following equation: Using the result from Lemma 3.1 and the Triangle Inequality, bearing in mind that ||e i || = 1 and that the basis field functions φ h are bounded, and after some algebraic manipulations, we can write the right-hand side as Coerciveness, see Lemma 2.2, and using φ h = w, gives hence there is an L 1 > 0 such that ||w|| H 1 (Ω) L (∆x) (n−1)/2 ε, which immediately implies that Theorem 3.2.Let u ε h be the solution to the boundary value problems in Eq (3.7), and u S h the solution to for all φ h ∈ V h (Ω). (3.9) Then there is an L 2 > 0 such that Proof.Using bilinearity of a(., .)gives upon setting w = u ε h − u S h the following equation: Using Taylor's Theorem for multivariate functions on smoothed delta distributions, we get the following result for the right-hand side: for x between x and x .The magnitude of the above expression can be estimated from above by We bear in mind that ∂ 3 δ With the two theorems above, we have proved that the solution to (BV P ε ) converges to the solution to (BV P ), and the solution to (BV P SP ) converges to the solution to (BV P ε ).Hence, we can derive the following theorem: Theorem 3.3.Let u h be the Galerkin solution to (BV P ) and u S h be the solution to (BV P SP ), let ε = O(∆x) p and ∆x → 0. If 0 < p < (2 + n)/3 then u S h converges to u h in the H 1 -norm, and u S h converges to u in the H 1 -norm.Proof.Denote u h and u S h to be the Galerkin solution to (BV P ε ) and (BV P SP ).Firstly, we consider From this inequality, we conclude that the finite element solution of the smooth particle method converges to the solution of the immersed boundary method upon letting ∆x → 0 and choosing ε = O(∆x) p ) for 0 < p < (2 + p)/3.

Numerical Results in Two Dimensions
To demonstrate the consistency between the immersed boundary approach and two alternative methods, we consider a square-shape cell in the computational domain.A homogeneous boundary condition is imposed for the exterior boundary of the computational domain.The parameter values are listed in Table 4.1.All of them are educated guesses in this study and they are dimensionless.
According to Proposition 3.1, to compare the immersed boundary approach and the 'hole' approach, the stiffness inside the biological cell needs to be adjusted, since two approaches   are consistent with β → 0. However, in the implementation, we can only select a very small positive value instead of β = 0. Numerical results are presented in Figure 4.1, Table 4.2 and Table 4.3.From the figure, there is no significant difference, except that in the smoothed particle approach, the displacement is a little bit larger than in the other two approaches.The reduction ratio of either the vicinity region or the cell appears to yield a tiny difference, which implies that three approaches are numerically consistent.However, the 'hole' approach takes slightly more computation time than the other two approaches.Therefore and due to the numerical complications in needing adaptive meshes, it will not be elected when we deal with the displacement and deformation of large number of cells, even though its convergence rate improves significantly comparing to the immersed boundary approach.As for the smoothed particle approach, the convergence rate of the L 2 −norm does not improve, while the computational efficiency does.

Conclusion
For the dimensionality exceeding one, the existence of Dirac Delta distributions in the elasticity equation results into a singular solution.We analyse the solutions based on Galerkin approximations with Lagrangian basis functions for different approaches that are consistent if cell sizes and smoothness parameters tend to zero.We have shown that all the alternative approaches are numerically consistent with the immersed boundary approach.The current paper has investigated and extended earlier findings to multi-dimensionality.The current analysis has been carried out for simple, linear elasticity.In the future, we plan to extend our findings to the visco-elasticity equations.This visco-elastic model contains a damping term, and still retains a linear nature.Furthermore, we are also interested in analyzing the above considered principles for a morphoelastic model.A morpho-elastic model has the major advantage of incorporating permanent deformations.A major complication is its nonlinear nature.

Figure 3 . 1 : 3 . 1 .
Figure 3.1: We consider a rectangular shape cell in two dimensions, with the centre position at (a, b).The forces exerted on the boundary are indicated by arrows (a) Immersed boundary approach (b) The 'hole' approach (c) Smoothed particle approach

Figure 4 . 1 :
Figure 4.1: For the different stiffness inside and outside of the cell, the solution (i.e. the displacement) is showed in each approach.Black curves show the deformed region of vicinity and the cell, and blue curve represents the cell.

Table 4 .
1: Parameter values used in the comparison of plasticity and morphoelasticity

Table 4 .
2:The percentage of area change of cell and vicinity region, and time cost of various approaches, if the stiffness is different inside and outside the biological cell.

Table 4 .
3: The L 2 − normof the solution (i.e. the displacement) with different mesh size in each approach, if the stiffness is different inside and outside the biological cell.