Upscaling between an agent-based model (smoothed particle approach) and a continuum-based model for skin contractions

Skin contraction is an important biophysical process that takes place during and after recovery of deep tissue injury. This process is mainly caused by fibroblasts (skin cells) and myofibroblasts (differentiated fibroblasts which exert larger pulling forces and produce larger amounts of collagen) that both exert pulling forces on the surrounding extracellular matrix (ECM). Modelling is done in multiple scales: agent-based modelling on the microscale and continuum-based modelling on the macroscale. In this manuscript we present some results from our study of the connection between these scales. For the one-dimensional case, we managed to rigorously establish the link between the two modelling approaches for both closed-form solutions and finite-element approximations. For the multi-dimensional case, we computationally evidence the connection between the agent-based and continuum-based modelling approaches.


Introduction
Wound healing is a spontaneous process for the skin to cure itself after an injury.It is a complicated combination of various cellular events that contribute to resurfacing, reconstituting and restoring of the tensile strength of the injured skin.For superficial wounds that only stay in epidermis, the wound can be healed without any problem.However, for severe injuries, in particular dermal wound, they may result in various pathological problems, such as contractures.
Contractures concur with disabilities and disfunctioning of patients, which cause a significantly severe impact on patients' daily life.Contractures are recognized as excessive and problematic contractions, which occur due to the pulling forces exerted by the (myo)fibroblasts on their direct environment, i.e. the extra cellular matrix (ECM).Contractions mainly start occurring in the proliferation phase of wound healing: the fibroblasts are migrating towards the wound and differentiating into myofibroblasts due to the high concentration of transforming growth factor-beta (TGF-beta).Usually, 5 − 10% reduction of wound area has been observed in clinical trials.A more detailed biological description can be found in [5,4,7,9].
In our previous work [13], a formalism to describe the mechanism of the displacement of the ECM has been used, which is firstly developed by Boon et al. [2] and improved further by Koppenol [8].Regarding the elasticity equation with point forces, we realized that the solution to the partial differential equation is singular for dimensionality exceeding one.Hence, we developed various alternatives to improve the accuracy of the solution in [11,12].
We have been working with agent-based models so far, which model the cells as individuals and define the formalism of pulling forces by superposition theory.However, once the wound scale is larger, the agent-based model is increasingly expensive from a computational perspective, and hence, the cell density model is preferred, which considers many cells as one collection in a unit.In this manuscript, we investigate and discover the connections between these two models, in the perspective of modelling the mechanism of pulling forces exerted by the (myo)fibroblasts.As the consistency between the smoothed particle approach (SPH approach) and the immersed boundary approach has been proven both analytically and numerically [11,12], we select the SPH approach here due to its continuity and smoothness, to compare with the cell density model using finite-element methods.
The manuscript is structured as follows.We start introducing both models in one dimension in Section 2, then in Section 3 we extend the models to two dimensions.Section 4 displays the numerical results in one and two dimensions.Finally, some conclusions are shown in Section 5.

Mathematical Models in One Dimension
Considering one-dimensional force equilibrium, the equations are given by Constitutive Equation.
By substituting E = 1, the equations above can be combined to Laplacian equation in one dimension:

Smoothed Particle Approach
In [12], a smoothed particle approach (SPH approach) is developed as an alternative of the Dirac Delta distribution describing the point forces exerted by the biological cells, in the application of wound healing: where P SP H is the magnitude of the forces, δ ε (x) is the Gaussian distribution with variance ε and s i is the centre position of biological cell i.One can solve the partial differential equations (PDEs) with finite-element methods.The corresponding weak form is given by Without this knowledge, the existence and uniqueness of the H 1 0 -solution follows as well from the application of the Lax-Milgram theorem [3], where it is immediately obvious that the bilinear form in the left-hand side is symmetric and positive definite.

Cell Density Approach
A cell density approach is often used in the large scale, so that the computational efficiency is much improved compared with the agent-based model.According to the model in [8], the force in two dimensions can be determined by the divergence of n c • I, where n c is the local density of the biological cells and I is the identity tensor.In one dimension, the cell density approach is expressed as: where P den is the magnitude of the forces.The corresponding weak form is given by for all φ ∈ H 1 0 ((0, L)).

Analytical Solutions with Specific Locations of Biological Cells
To express the analytical solution, it is necessary to determine the locations of the biological cells, such that the cell density can be written as an analytical function of the positions.We assume, there are N s cells distributed uniformly in the subdomain (a, b) of the computational domain (0, L).Hence, the distance between the center position of any two adjacent biological cells is constant, which we denote ∆s = (b − a)/N s and the first and the N s -th cell are located at x = a + ∆s/2 and x = b − ∆s/2, respectively.With homogeneous Dirichlet boundary conditions, and suppose P SP H = P ∆s and variance ε = ∆s, the boundary value problem of the SPH approach is expressed as where P is a positive constant and s i is the centre position of the biological cells.Utilizing the superposition principle, the analytical solution is given by where erf(x) is the error function defined as erf(x) = 2 √ π x 0 exp(−t 2 )dt [15].Since the biological cells are uniformly located between a and b (0 < a < b < L), dnc dx can be rephrased as where t is a small positive constant.Taking t to zero, the above expression converges to δ(x − a) − δ(x − b).Hence, the boundary value problem of the cell density model can be written as where δ(x) is the Dirac Delta distribution and a and b are the left and right endpoint of the subdomain (where biological cells are uniformly located) respectively.The analytical solution is then expressed as where G(x, x ) is the Green's function [6], defined by in the computational domain (0, L).
Actually, the convergence between u 1 (x) and u 2 (x) can be proven as ∆s → 0 + by a proposition.Firstly, we introduce Chebyshev's Inequality: Lemma 2.1.(Chebyshev's Inequality [10]) Denote X as a random variable with finite mean µ and finite variance σ 2 .Then for any positive k ∈ R, the following inequality holds: where P(A) is the probability of event A. The above inequality can also be rephrased as Proposition 2.1.Let u 1 (x) as described in Eq (2.4) be the exact solution to (BV P 1 SP H ) and u 2 (x) as described in Eq (2.6) be the exact solution to (BV P 1 den ).As ∆s → 0 + , u 1 (x) converges to u 2 (x).
Proof.For the standard Gaussian distribution in one dimension, the cumulative distribution function is given by F . Thus, we obtain (2.8) By Chebyshev's Inequality (see Lemma 2.1), one can conclude that for any positive k, (2.9) Note that 1 − F (k) = F (−k) due to the symmetry of standard Gaussian distribution.Hence, Eq (2.9) implies and analogously, 0 As it has been defined earlier that ∆s = (b − a)/N s , we obtain Analogously, we obtain that for any Thus, it can be concluded that for any series of real number {x i } ∈ R n , when x i+1 − x i = ∆s and x i is either all positive or all negative for any i where sgn(x) is sign function defined by We rewrite u 1 (x) as Combining Eq (2.10), (2.11) and (2.12), u 1 (x) is given by Rewriting u 2 (x) regarding different domain gives Hence, we conclude that u 1 (x) converges to u 2 (x) as ∆s → 0 + .

Finite-Element Method Solutions with Arbitrary Locations of Biological Cells
For the finite-element method, we select the piecewise Lagrangian linear basis functions.We divide the computational domain into N e mesh elements, with the nodal point x 1 = 0 and x Ne+1 = L.For the implementation, we define the cell density as the count of biological cell in every mesh element divided by the length of the mesh element, hence, it is a constant within every mesh element.In other words, in the mesh element [x j , x j+1 ], the count of the biological cell is defined by for any j ∈ {1, . . ., N e }, where h is the size of every mesh element.Different from (BV P 1 SP H ) where ∆s is the variance of δ ε , for finite-element methods, we set ε = h/3, such that the integration of δ h/3 (x − x ) for any 0 < x < L over any mesh element with size h, is close to 1 (see Lemma 2.2).With the two approaches, the boundary value problems with Dirichlet boundary condition are defined by and where s i is the position of biological cells, h is the mesh size and N s is the total number of cells in the computational domain.The consistency between (BV P 2 SP H ) and (BV P 2 den ) can be verified by the following lemma and theorem.Lemma 2.2.(Empirical rule [14]) Given the Gaussian distribution of mean µ and variance ε: then the following integration can be computed: respectively the solution to (BV P 2 SP H ) and (BV P 2 den ).With Lagrangian linear basis functions for the finite element method, u h 1 (x) converges to u h 2 (x), as the size of the mesh element h → 0 + , regardless of the positions of biological cells.
Proof.We define v h (x) = u h 1 (x) − u h 2 (x), then the boundary value problem to solve v h (x) is given by (2.15) The corresponding Galerkin's form reads as for all φ ∈ H 1 0 ((0, L)).Using integration by parts and letting φ = φ j , j ∈ {1, . . ., N e + 1}, the equation in (GF 1 v ) can be rewritten by 3 Mathematical Models in Two Dimensions

Smoothed Particle Approach and Cell Density Approach
In multi dimensional case, the equation of conservation of momentum over the computational domain Ω, without considering inertia, is given by We consider a linear, homogeneous and isotropic domain, with Hooke's Law, the stress tensor σ is defined as where E is the Young's modulus of the material, ν is Poisson's ratio and is the infinitesimal strain tensor: Considering a subdomain Ω w ⊂ Ω, where the center positions of the biological cells are located, then the SPH approach and cell density approach with homogeneous Dirichlet boundary condition are derived by and

Consistency between Two Approaches in Finite-Element Method
To prove the consistency between these two approaches, we define that for the triangular mesh element e k , k ∈ {1, . . ., N e }, where N e is the total number of mesh elements in Ω, the density of biological cell n c (e k ) is constant and the count of biological cells N c (e k ) is expressed by where A(e k ) is the area of mesh element e k .Similar to the process of proof in one dimension, we state the following theorem: Theorem 3.1.Denote u h 1 (x) and u h 2 (x) respectively the solution to (BV P 3 SP H ) and (BV P 3 den ) with P SP H = P den = P .With Lagrange linear basis functions for the finite element method, u h 1 (x) converges to u h 2 (x), as the size of a triangular mesh element is taken to zero, regardless of the positions of biological cells.
Proof.We consider v h (x) = u h 1 (x) − u h 2 (x), then the boundary value problem to solve v h (x) is given by The corresponding Galerkin's form reads as With integral by parts and letting φ h = φ h k , k ∈ {1, . . ., N e }, the equation in (GF 3 v ) can be rewritten by δ(x − s i )dΩ.Note that in two dimensions, ε needs to be sufficiently small compared to the size of a triangular mesh element.

Results
Results in both one and two dimensions are discussed in this section.Since the objective of this manuscript is to investigate the consistency and the connections between the SPH approach and the cell density approach, all the parameters are dimensionless.

One-Dimensional Results
We show the results by analytical solutions in Figure 4.1 with various values of ∆s (i.e.depending on different number of biological cells in the subdomain (a, b)).Here, the computational domain is (0, 7) with L = 7 and the subdomain where the biological cells locate uniformly is (2, 5) with a = 2 and b = 5.With the decrease of the variance in the Gaussian distribution in (BV P 2 SP H ), the curves gradually overlap, which verifies the convergence between the analytical solutions to these two approaches.To implement the model, there are two different algorithms shown in Figure 4.2 and 4.3.Depending on different circumstances, the implementation method is elected.The cell density in one dimension is defined as the number of cells per length unit.In other words, the cell count in a given domain can be computed by integrating the cell density over the domain.If the cell density function can be expressed analytically and the first order derivative of the function exists, then a certain bin length d is chosen and the cell count in every bin of d length is calculated.Then we generalize the center positions of cells in every bin of length d, thus, the SPH approach can be implemented, as it is indicated in Figure 4.2.However, it is not always straightforward to obtain the analytical expression of cell density.If the center positions of cells are given, the number of cells in each mesh element can be counted, hence, the cell density will be computed analogously at each mesh points.Therefore, the boundary value problem of cell density approach is solved by numerical methods, for example, the finite-element methods.

Cell density approach
Cell positions SPH approach Figure 4.2: With exact expression of cell density function and the first order derivative of the function exists, cell density approach is implemented directly.Based on the cell density, the number of cells in a certain region with length d is determined and subsequently, the center positions of cells can be generalized.Hence, the SPH approach is implemented.

Cell positions
Cell count in mesh element

SPH approach
Cell density Cell density approach Computing the number of cells in every mesh element and divided by the length of the mesh element results into the cell density.Subsequently, cell density approach can be implemented.
In this manuscript, all the numerical results are derived by finite-element methods with Lagrangian linear basis functions.Regarding the first implementation method (see Figure 4.2), we show the results with a Gaussian distribution and sine function as cell density functions; see Figure 4.4 and 4.5.We start with the simulations in which we keep the number of cells and the center positions of the cells the same, then we refine the mesh.In Figure 4.5(a)-(c), the bin length d is 0.35, and the mesh size is a function of d.The results solved by SPH approach become smoother.With various values of d, the solutions to the approaches are overlapping only when the factor between the d and mesh size is closer to 1. From Figure 4.5(d) to (f), the mesh is fixed and we vary the value of d.We note that in Figure 4.5(f), the solution to the SPH approach is significantly different from the solution to the cell density approach.It is mainly caused by the fact that d is too small and there is barely any fluctuation with the count of cells in every d length subdomain, while with the Gaussian distribution as the cell density function, the majority of the cells are centered around x = 3.5.Hence, the solution to SPH approach still manages to be comparable with the solution to the cell density approach; see Figure 4.4(f).Numerical results of the simulation in Figure 4.4 are displayed in Table 4.1.There are some noticeable differences between two approaches, in particular the convergence rate in the H 1 -norm: thanks to the given, differentiable cell density function, the cell density approach converges faster.In addition, the cell density approach requires less computational time with a factor of 15.
Here, we define N s = 88 and the mesh size h = 0.07.The results are solved by finite-element method with algorithm in Figure 4.2.

SPH Approach
Cell Density Approach We consider cells that are located uniformly in the subdomain (2,5), which implies that the gradient or divergence of the cell density vanishes inside the subdomain but does not exist at two endpoints of the subdomain.Hence, we utilize the implementation method in    show the solutions to (BV P 2 SP H ) and (BV P 2 den ) respectively.Note that, in the finite-element method solutions, the magnitude of the forces in both approaches are the same, and the variance of δ ε (x) is related to h rather than ∆s.Furthermore, these figures verify that the convergence between SPH approach and cell density approach is determined by the mesh size rather than by the distance between any two adjacent cells.Table 4.2 displays the numerical results of the simulation in Figure 4.6, in the perspective of the solution, the reduction ratio of the subdomain and the computational cost.Similarly to the figures, there is no significant difference between the norms and the deformed length of the subdomain.However, the simulation time in the cell density approach is much shorter than in the SPH approach with a factor of 35.

Two-Dimensional Results
In the multi-dimensional case, we are not able to write the analytical solution to the boundary value problems.The results are all solved by the use of the finite-element method applied to (BV P 3 SP H ) and (BV P 3 den ).Note that the force magnitude of both boundary value problems is the same.Following the same implementation methods as in one dimension, simulations are carried out with two formulas of cell density: (1)  For Simulation (1), no analytical expression for (the derivative of) the density function is available.Therefore, the implementation starts with generating the cell positions, according to the principles outlined in Figure 4.3.In Figure 4.8, the displacement results are displayed.From the figures, hardly any significant differences between the solutions can be observed, which indicates that these two approaches are numerically consistent.Table 4.3 displays more details about the two approaches regarding the numerical analysis: most data are more or less the same.Thanks to the continuity of the SPH approach, the convergence rate between two approaches is similar.However, as it has been mentioned earlier, the agent-based model is computationally more expensive than the continuum-based model; here, the difference is a factor of 240.
According to the setting of the simulation, we define the cell density function by  which is a Gaussian distribution multiplied by a positive constant.Similarly, in two dimensions, the cell density is defined by the number of biological cells per unit area.In other words, the cell count is computed by the local cell density multiplied by the area of selected region.Here, we assume that the selected region is a 1 × 1 square, then we generate the center positions of biological cells in every unit square based on the local number of cells.Figure 4.9 shows the numerical results regarding two approaches.There is no significant difference if the same mesh structure resolution is used.As the mesh is refined, the solution to the SPH approach is smoother, since the "ring" in the center becomes more dominant.In Table 4.4, it can be concluded that there are no significant differences between two approaches except for the computational efficiency and the convergence rate of H 1 -norm.If the mesh is not fine enough, then the solution to the SPH approach is less smooth, hence, the determination of the gradient of the solution is less accurate.

Conclusions
In this manuscript, we discussed the different models to simulate the pulling forces exerted by the (myo)fibroblasts depending on different scales of the wound region.We started from one dimension and later extended the models to two dimensions.In one dimension, we can write explicitly the analytical solution to the boundary value problem with specific distribution of the locations of biological cells.The consistency of the solutions indicates the possibility to prove the analytical connection between the two approaches.In both one and two dimensions, the numerical solutions delivered by finite-element methods with Lagrange linear basis functions implied that these two models are consistent under certain mesh conditions (when the mesh size is sufficiently small) and regardless the locations of the biological cells and the implementation methods.In summary, regarding the displacement of the ECM from the mechanical model, the agent-based model and the cell density model are consistent from a computational point of view.This could be used to transfer one type of model to the other one regarding the force balance in the wound healing model, as the connection between these two models has been suggested.We want to use the developed insights for the analysis of upscaling between agent-based and continuum-based model formulations.

Figure 4 . 1 :
Figure 4.1: The exact solutions to (BV P 1 SP H ) and (BV P 1 den ) are shown, with various values of ∆s, which is the distance between centre positions of any two adjacent biological cells.Blue points are the centre positions of biological cells.Red curves represent the solutions to (BV P 1 SP H ) and blue curves represent the solutions to (BV P 1 den ).

Figure 4 . 3 :
Figure 4.3: Given the center positions of cells, one can directly implement the SPH model.Computing the number of cells in every mesh element and divided by the length of the mesh element results into the cell density.Subsequently, cell density approach can be implemented.
Figure 4.3, as the center positions of the cells are given, then the local cell density can be calculated per unit area.Compared with the results shown in Figure 4.1, the results in Figure 4.6 and Figure 4.7

0175 Figure 4 . 4 :
Figure 4.4: The cell density function is Gaussian distribution and using the algorithm in Figure 4.2, different simulations are carried out with various mesh size and the total number of cells.Blue curves represent the solutions to (BV P 2 SP H ), and ref curves are the solutions to (BV P 2 den ) with n c (x) = 50 × 1/ √ 2π × 0.1 2 exp{−(x − 3.5) 2 /(2 × 0.1 2 )}.In Subfigure (a)-(c), we set d = 0.35 and cell positions are fixed, as h is decreasing.From Subfigure (d) to Subfigure (f ), we use the same finite-element method settings (where h is sufficiently small with h = 0.07), and simulations are carried out with various values of d.

Figure 4 . 5 :
Figure 4.5: The cell density function is sine function and using the algorithm in Figure 4.2, different simulations are carried out with various mesh size and the total number of cells.Blue curves represent the solutions to (BV P 2 SP H ), and ref curves are the solutions to (BV P 2 den ) with n c (x) = 40| sin(2x)|.In Subfigure (a)-(c), we set d = 0.35 and cell positions are fixed.From Subfigure (d) to Subfigure (f ), we use the same finite-element method settings (where h is efficiently small with h = 0.07), and we take different values of d.

Figure 4 . 6 :
Figure 4.6: The finite-element method solutions to (BV P 2 SP H ) and (BV P 2 den ) are shown where cells are uniformly located.With the fixed positions of cells, the solutions are convergent as h → 0 + .Blue points are the centre positions of biological cells.Red curves represent the solutions to (BV P 2 SP H ) and blue curves represent the solutions to (BV P 2 den ).

Figure 4 . 7 :
Figure 4.7: The finite-element method solutions to (BV P 2 SP H ) and (BV P 2 den ) are shown with uniform distribution.Compared to the analytical result, the consistency between two approaches are unrelated to the number of cells, and the solutions are convergent as h → 0 + .Here, we use h = 0.007.Blue points are the centre positions of biological cells.Red curves represent the solutions to (BV P 2 SP H ) and blue curves represent the solutions to (BV P 2 den ).

2 .
cells are located inside the subdomain Ω w randomly by the uniform distribution; (2) the cell density function is in the form of the standard Gaussian distribution over the computational domain with n c (x) = 50 × 1 2π exp − x| 2 Implementation methods in Figure 4.3 and Figure 4.2 are applied respectively in Simulation (1) and (2).

Figure 4 . 8 :
Figure 4.8: The displacement results are shown, which are solved from (BV P 3 SP H ) and (BV P 3 den ) when cells are randomly located in the subdomain (−5, 5) × (−5, 5).In other words, it is impossible to write the analytical expression of n c (x), subsequently, the algorithm in Figure 4.3 is selected.There are 196 biological cells in the computational domain.Blue points are the center positions of biological cells, red curves are the original shapes of the subdomain, and the black curves represent the deformed boundary of the subdomain.

Figure 4 . 9 :
Figure 4.9: The displacement results are shown, which are solved from (BV P 3 SP H ) and (BV P 3 den ) when cells are located according to the cell density function n c (x) = 50 × 1 2π exp{− x| 2 2 }.Hence, implementing algorithm in Figure 4.3 is used.There are 440 biological cells in the computational domain.Blue points are the center positions of biological cells, red curves are the original shapes of the subdomain, and the black curves represent the deformed boundary of the subdomain.

Table 4 .
1: Numerical results of two approaches in one dimension, where the cell density function is Gaussian distribution: n

Table 4 .
2: Numerical results of two approaches in one dimension with biological cells located uniformly.Here, we define mesh size h = 0.07 and N s = 50, which means ∆s = 0.06.The results are solved by finite-element method with algorithm in Figure4.3.

Table 4 . 3 :
Numerical results of two approaches in two dimensions with random distribution for the positions of biological cells.Due to the nonexistence of divergence or gradient of cell density function, implementation method in Figure4.3 is used.

Table 4 .
4: Numerical results of two approaches in two dimensions with Gaussian distribution for the positions of biological cells.Figure4.2 is implemented and there are 440 biological cells in the computational domain.