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 processes that contribute to resurfacing, reconstituting and restoring of the tensile strength of the injured skin. For superficial wounds where only the epidermis is damaged, they heal without any issues. However, for severe injuries, in particular dermal wounds, they may result in various pathological problems, such as contractures.

Contractures concur with disabilities and disfunctioning of joints 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, reduction of wound area has been observed in clinical trials. A more detailed biological description can be found in Enoch and Leaper (2008); Cumming et al. (2009); Haertel et al. (2014); Martin (1997).

In our previous work (Peng and Vermolen 2020b), a formalism to describe the mechanism of the displacement of the ECM has been used, which is firstly developed by Boon et al. (2016) and improved further by Koppenol (2017). To model skin contraction, which results from the cellular traction forces applied on the ECM, the momentum balance equation (i.e. the elasticity equation) is combined with the immersed boundary approach (Ferziger et al. 2002). Regarding the elasticity equation with point forces, we realized that the solution to the partial differential equation is singular in the sense that the formal solution does not reside in the same function space as the finite element solution does for two- and higher dimensional problems. Hence, we developed various alternatives to improve the accuracy of the solution in Peng and Vermolen (2020c, 2022).

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, which is a continuum-based model, is preferred. 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. Since the consistency between the smoothed particle approach (SP approach) and the immersed boundary approach has been proven both analytically and numerically (Peng and Vermolen 2019a, b), we select the SP approach here due to its continuity and smoothness, to compare with the cell density model using finite-element methods.

The objective of this paper is to demonstrate the consistency between the smoothed particle approach and the cell density approach, therefore, the same physical scale is considered in the current study. The cell density approach, which amounts to the continuum scale is mainly considered in scales that are in the order millimeters or larger. For such scales, the number of cells becomes too large and hence it is no longer computationally feasible to use individual cells that have sizes that range in the order of micrometers or centimeters and hence a continuum-based approach with cell densities is necessary so that cellular processes such as migration (random walk, chemotaxis etc.), cell division/proliferation, cell death can be dealt with by the incorporation of several terms in partial differential equations. For small scales, ranging in magnitudes that are smaller than millimeters, averaged cell densities no longer make sense and hence the continuum-based models become less useful and hence partial differential equations can no longer approximate the dynamics of cellular processes. Another reason for the use of agent-based models is their formulation in terms of measurable quantities such as cell forces, cell migration velocities etc. For even smaller scales that range within nanometers, the partial differential equations that we employ for the balance of momentum is no longer applicable. Then one has to use molecular dynamics-based simulations for all processes, including the mechanics.

The manuscript is structured as follows. We start introducing both models in one dimension in Sect. 2, then in Sect. 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 Sect. 5.

Mathematical models in one dimension

This section considers one-dimensional solutions for two models: the smoothed particle model and the cell density model. Furthermore, convergence of the analytic solutions of these approaches, as well as convergence of the numerical solution are considered in this chapter. The treatment of the one-dimensional case is relatively straightforward due to the simple nature of the exact solutions.

The balance of momentum, is modelled by linear elasticity based on isotropy. In the Navier-Cauchy equation, inertia is neglected. For a general dimensionality, the Navier-Cauchy equation consists of Eqs. (16), (17) and (18) that are displayed in Sect. 3. We consider the one-dimensional version of the Navier-Cauchy equation in an isotropic and continuous domain, hence, the equations are given byBy substituting , the equations above can be combined to Poisson equation in one dimension:

Smoothed particle approach

In Peng and Vermolen (2019b), a smoothed particle (SP) 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. By specifying the force expression f in Eq. (1) and considering cells, the smoothed particle approach (Peng and Vermolen 2020a, c, 2021) is given bywhere is the magnitude of the forces, is the Gaussian distribution with variance and is the centre position of biological cell i. The boundary conditions close the problem so that it admits a uniquely defined solution. One can solve the partial differential equations (PDEs) with finite-element methods. The corresponding weak form is given byThe existence and uniqueness of the -solution follows as well from the application of the Lax–Milgram theorem (Braess 2007), 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 Koppenol (2017), the force in two dimensions is proportional to the divergence of , where is the local density of the biological cells and is the identity tensor. In one dimension, this becomes , the cell density approach is expressed as:where is the magnitude of the forces. The corresponding weak form is given by

Consistency between two models

In the application of wound healing, we assume an artificial wound embedded within the computational domain. Therefore, for the one-dimensional case, we define the computational domain as (0, L) as aforementioned, and biological cells are located in the subdomain where ; see Fig. 1 for a schematic representation.

Fig. 1

A schematic representation of the computational domain (0, L) in one dimension. The subdomain with where the biological cells (red dots in the figure) are only located is marked as a blue bar

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 that there are cells distributed uniformly in the subdomain (a, b) of the computational domain (0, L). Hence, the distance between the center positions of any two adjacent biological cells is constant, which we denote by and the first and the -th cell are located at and , respectively. With homogeneous Dirichlet boundary conditions, and suppose and variance , the boundary value problem of the SP approach is expressed aswhere P is a positive constant and is the centre position of the biological cells. Utilizing the superposition principle, the analytical solution (i.e. the displacement at arbitrary position of the domain) is given bywhere is the error function defined as (Weisstein 2010). Note that the solution satisfies the Dirichlet boundary conditions in .

Since the biological cells are uniformly located between a and b (), can be rephrased aswhere t is a small positive constant. In other words, increases linearly in and decreases linearly in — with respect to x — and stays constant elsewhere 1/t in and zero for . Taking t to zero, the above expression converges to . Hence, the boundary value problem of the cell density model can be written aswhere 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 (i.e. the displacement at arbitrary position of the domain) is then expressed aswhere is the Green’s function (Haberman 1983), defined byin the computational domain (0, L). Note that the solution satisfies the Dirichlet boundary conditions in .

We will demonstrate the convergence between and as . First, we introduce Chebyshev’s Inequality:

Lemma 1

(Chebyshev’s Inequality (Olkin and Pratt 1958)) Denote X as a random variable with finite mean and finite variance . Then for any positive , the following inequality holds:where is the probability of event A. The above inequality can also be rephrased as

The proof of Chebyshev’s Inequality is standard (Olkin and Pratt 1958), and therefore we do not give it here.

Proposition 1

Let as described in Eq. (4) be the exact solution to and as described in Eq. (6) be the exact solution to . As , converges to .

Proof

For the standard Gaussian distribution in one dimension, the cumulative distribution function is given by . Thus, we obtainBy Chebyshev’s Inequality (see Lemma 1), one can conclude that for any positive k,Note that due to the symmetry of standard Gaussian distribution. Hence, Eq. (9) impliesand analogously, is implied. Together with Eq. (8), it givesLet , for any , thenAs it has been defined earlier that , we obtainSince for any , the Squeeze Theorem (Apostol and Ablow 1958) implies thatAnalogously, we obtain that for any ,Thus, it can be concluded that for any series of real number , when and is either all positive or all negative for any ,where is the sign function defined byWe rewrite asCombining Eqs. (10), (11) and (12), is given byRewriting regarding a general bounded domain givesHence, we conclude that converges to as .

Finite-element method solutions with arbitrary locations of biological cells

For the finite-element method, we use piecewise Lagrangian linear basis functions. We divide the computational domain into mesh elements, with the nodal point and . For the implementation, we define the cell density as the count of biological cells in every mesh element divided by the length of the mesh element, hence, it is a constant within every mesh element and it is an interval function over the mesh elements. In other words, in the mesh element , the count of the biological cell is defined byfor any , where h is the size of every mesh element. To clarify the notations, we use for the cell density function when it is an interval function in one dimension, and when it can be written analytically as a continuous function. Same settings hold for the function of cell count . Different from where is the variance of , for finite-element methods, we set , such that the integration of for any over any mesh element with size h, is close to 1 (see Lemma 3, which follows later). With the two approaches, the boundary value problems with Dirichlet boundary condition are defined byandwhere is the position of biological cells, and is the total number of cells in the computational domain. The consistency between and can be verified by the following lemmas and theorem.

Lemma 2

Given the Gaussian distribution of mean and variance :then for any , it follows that

We use the transformation , which transforms the bounds of the integral, and further changes the integral intoSending , makes the bounds tend to infinity, hencesince it is well known thatThis concludes the proof.

In order to quantify the error of having a value differing from zero, one can use the following Empirical Rule:

Lemma 3

(Empirical rule (Pukelsheim (1994))) Given the Gaussian distribution of mean and variance :then the following integration can be computed:

Theorem 1

Denote and respectively the finite-element solution to and . With Lagrangian linear basis functions for the finite element method, converges to , as in the Gaussian distribution of the SP approach, regardless of the positions of biological cells.

We define , then satisfies the following boundary value problemThe corresponding Galerkin’s form reads asUsing integration by parts and letting , the equation in can be rewritten bysince the number of biological cells in for any can be defined as , where is the total number of biological cells in the computational domain. We note that if a cell center is located on a nodal point, such as or , then only half of the unit counts as . We, further, note that Empirical Rule, Lemma 3, can be used to quantify the error for not being identical zero.

Mathematical models in two dimensions

Smoothed particle approach and cell density approach

In the multi dimensional case, the Navier-Cauchy equation of conservation of momentum over the computational domain which is an isotropic and homogeneous domain, without considering inertia, is given bySince we consider a linear, homogeneous and isotropic domain, with Hooke’s Law, the stress tensor is defined aswhere E is the Young’s modulus of the material, is Poisson’s ratio and is the infinitesimal strain tensor:Considering a subdomain , where the center positions of the biological cells are located, then the SP approach and cell density approach with homogeneous Dirichlet boundary condition are derived byandwhere is the unknown displacement field which is to be found.

Consistency between two approaches in finite-element method

To prove the consistency between these two approaches, we define that for the triangular mesh element , where is the total number of mesh elements in , the density of biological cells is constant within a mesh element and a function of the mesh element (hence it is a piecewise constant function), and hence the count of biological cells is expressed bywhere is the area of mesh element . To indicate that the Gaussian distribution is a proper replacement of the Dirac Delta distribution as in two dimensions, we state the following lemma:

Lemma 4

Given the bivariate Gaussian distribution centered at mean with variance :then for any ball centered at with radius R, denoted by , it follows thatHence, for any , it follows that

Proof

We use centered polar coordinates, and , and the Jacobian of the transformation, to arrive atTreating as an arbitrary constant, and sending to zero, givesThis concludes the proof.

Then, similar to the one dimensional case, we state the following theorem:

Theorem 2

Let and , respectively, denote the finite-element solution to and with . With Lagrangian linear basis functions for the finite element method, converges to , as in the Gaussian distribution of the SP approach, regardless of the positions of biological cells.

We consider , satisfies the following boundary value problem:The corresponding Galerkin’s form reads asWith integration by parts and letting , the equation in can be rewritten bysince it can be defined that , where is the total number of biological cells in the computational domain. Note that in two dimensions, needs to be sufficiently small compared to the size of a triangular mesh element. Korn’s inequality (Braess 2007) and symmetry (boundedness) conclude the theorem. We note that if cell is located on a face of an element, then its contribution only counts for a half as . If a cell is located at an element vertex, then its contribution only counts for the angle of the vertex in the current element divided by as .

We note that Lemma 4 can be used to estimate the error when one wishes to accommodate for not being identically zero.

Simulation results

Simulation 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 SP approach and the cell density approach, all the parameters are dimensionless. For one and two dimensional simulations, the parameter values are shown in Tables 1 and 5, respectively.

Table 1

Model parameter values that are used in one dimensional simulations

Parameter	Value	Description
E	1	Stiffness of the computational domain (i.e. ECM)
P	0.01	Force magnitude in both approaches as indicated in (BVPv1) in Theorem 1
L	7	The right endpoint of the computational domain
a	2	The left endpoint of the subdomain (i.e. the wound region)
b	5	The right endpoint of the subdomain (i.e. the wound region)

Table 5

Estimated parameter values that are used in one dimensional simulations

Parameter	Value	Description
E	1	Stiffness of the computational domain (i.e. ECM)
P	10	Force magnitude in both approaches as indicated in (BVPv3) in Theorem 2
ν	0.49	Poisson’s ratio of the computational domain (i.e. ECM)
x_0	20	The length of the computational domain in x-coordinate
y_0	20	The length of the computational domain in y-coordinate
wx	10	The length of the subdomain (i.e. the wound region) in x-coordinate
wy	10	The length of the subdomain (i.e. the wound region) in y-coordinate

One-dimensional results

We show the results by analytical solutions (see Eqs. (5) and (7) for the SP approach and the cell density approach, respectively) in Fig. 2 with various values of (i.e. depending on different number of biological cells in the subdomain (a, b)). Here, the computational domain is (0, 7) with and the subdomain where the biological cells locate uniformly is (2, 5) with and . With the decrease of the variance in the Gaussian distribution in , the curves gradually overlap, which verifies the convergence between the analytical solutions to these two approaches.

Fig. 2

The exact solutions to and are shown, with various values of , 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 and blue curves represent the solutions to

To implement the model, there are two different algorithms shown in Figs. 3 and 4. 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 SP approach can be implemented, as it is indicated in Fig. 3. However, it is not always straightforward to obtain the analytical expression of cell density. If the center positions of cells are given, then the SP approach can be implemented directly, and the number of cells in each mesh element can be counted. Hence, the cell density will be computed analogously at each mesh points, as it is shown in Fig. 4. Therefore, the boundary value problem of cell density approach is solved by numerical methods, for example, the finite-element methods. In summary, in Fig. 3 when the cell density is prescribed and on the basis of this cell density function, we assign a number of cells in each line element (in one dimension) or triangular (in two dimensions) element. This amounts to using a constant cell density in each element (line element or triangular element). The approach shown in Fig. 4 is the contrary, where cell density cannot be expressed analytically and then the number of cells in each element is counted and divided by the correspond measure (length in one dimension and element area in two dimensions) to reconstruct the cell density distribution. In both cases, the cell density is piecewise constant over the domain of computation in which the cell density is constant over each element.

Fig. 3

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 SP approach is implemented

Fig. 4

Given the center positions of cells, one can directly implement the SP 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

In this manuscript, all the numerical results have been obtained by finite-element methods with Lagrangian linear basis functions. Regarding the first implementation method (see Fig. 3), we show the results with a sine function and Gaussian distribution as cell density functions; see Figs. 5 and 6 , respectively. 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 Fig. 5(a)–(c), the bin length d is 0.35, and the mesh size is a function of d. The results that were obtained using the SP 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 Fig. 5(d) to (f), the mesh is fixed and we vary the value of d. We note that in Fig. 5(f), the solution to the SP approach is significantly different from the solution to the cell density approach. This difference is mainly caused by the fact that d is too small and there is barely any fluctuation with the count of cells in every subdomain with length d, while with the Gaussian distribution as the cell density function, the majority of the cells are centered around . Hence, the solution that was obtained from the SP approach still manages to be comparable with the solution to the cell density approach; see Fig. 6(f). Numerical results of the simulation in Fig. 6 are displayed in Table 2. There are some noticeable differences between two approaches, in particular the convergence rate in the -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.

Fig. 5

The cell density function is sine function and using the algorithm in Fig. 3, different simulations are carried out with various mesh size and the total number of cells. Blue curves represent the solutions to , and red curves are the solutions to with . In (a)–(c), we set and cell positions are fixed. From (d) to (f), we use the same finite-element method settings (where h is efficiently small with ), and we take different values of d

Fig. 6

The cell density function is Gaussian distribution and using the algorithm in Fig. 3, different simulations are carried out with various mesh size and the total number of cells. Blue curves represent the solutions to , and red curves are the solutions to with . In (a)–(c), we set and cell positions are fixed, as h is decreasing. From (d) to (f), we use the same finite-element method settings (where h is sufficiently small with ), and simulations are carried out with various values of d

Table 2

Numerical results of two approaches in one dimension, where the cell density function is Gaussian distribution: . Here, we define and the mesh size . The results are solved by finite-element method with algorithm in Fig. 3

	SP Approach	Cell Density Approach
L^2-norm of the solution u
‖u‖L^2((0,L))2h	0.5419891028287082	0.3581360438718032
‖u‖L^2((0,L))h	0.5441481069175041	0.36197930815501245
‖u‖L^2((0,L))h/2	0.5448668763153627	0.3630629995429662
Convergence rate of L^2-norm	1.75281178	1.826378221
‖u‖H^1(((0,L))h	0.9642151731656272	0.871720645462775
Reduction ratio of the subdomain (a, b) (%)	13.88062	9.52381
Relative ratio of the subdomain (a, b) after deformation (%)	86.11938	90.47619
Time cost (s)	0.045070	0.0032084

We consider cells that are located uniformly in the subdomain (2, 5), which implies that the derivative of the cell density vanishes inside the subdomain but does not exist at two endpoints of the subdomain. Since the exact solutions of both approaches are known as in Eq. (5) and in Eq. (7), respectively, one can perform root-mean-square error (RMS) analysis, which is given bywhere N is the number of mesh nodal points and is the coordinate of th mesh nodal point. Note that the error above is actually error computed in the computational domain. To obtain the numerical results and since the cell density can be written analytically, we utilize the implementation method in Fig. 4, 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 Fig. 2, the results in Figs. 7 and 8 show the solutions to and respectively. Note that, in the finite-element method solutions, the magnitude of the forces in both approaches are the same, and the variance of is related to h rather than . Furthermore, these figures verify that the convergence between SP approach and cell density approach is determined by the mesh size rather than by the distance between any two adjacent cells. Table 3 displays the numerical results of the simulation in Fig. 7, 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 SP approach with a factor of 35.

Fig. 7

The finite-element method solutions to and are shown where cells are uniformly located. With the fixed positions of cells, the solutions are convergent as . Blue points are the centre positions of biological cells. Red curves represent the solutions to and blue curves represent the solutions to

Fig. 8

The finite-element method solutions to and 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 . Here, we use . Blue points are the centre positions of biological cells. Red curves represent the solutions to and blue curves represent the solutions to

Table 3

Numerical results of two approaches in one dimension with biological cells located uniformly. Here, we define mesh size and , which means . The results are solved by finite-element method with algorithm in Fig. 4

	SP Approach	Cell Density Approach
RMS error	6.70984672×10^-7	0.010295901
L^2-norm of the solution u
‖u^h‖_L^2((0,L))	0.21470072397236404	0.2190152169099139
‖u^h/2‖_L^2((0,L))	0.2180816688838546	0.2192106132836521
‖u^h/4‖_L^2((0,L))	0.21897037944459403	0.219296655339645
Convergence rate of L^2-norm	1.927640972	1.831863826
Reduction ratio of the subdomain (a, b) (%)	7.96908	7.98821
Relative ratio of the subdomain (a, b) after deformation (%)	92.03092	92.01179
Time cost (s)	0.10391	0.0030458

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 and . 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) the cell density function is in the form of the standard Gaussian distribution over the computational domain with ; (2) cells are located inside the subdomain randomly by the uniform distribution. Implementation methods in Figs. 3 and 4 are applied respectively in Simulation (1) and (2).

According to the setting of the simulation, we define the cell density function bywhich 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 square, then we generate the center positions of biological cells in every unit square based on the local number of cells.

In the two-dimensional calculations (see Fig. 9), we consider a square domain for the tissue. Within this domain, there is an artificial scar, which is also square-shaped in the current simulations, and which is indicated by the red line segments. The scar, which is bounded by the red line segments, is populated with cells (fibroblasts) that exert pulling forces. The cells are positioned randomly in the scar region such that they do not overlap. In the smoothed particle approach, the gradient of the mollified Dirac delta distribution, which amounts to a Gaussian distribution, is used to model the force exerted by each cell, whereas in the cell density approach a cell density field is reconstructed from the randomized cell positions using the procedure outlined in and ; see Sect. 3.1. Summarized, in both cases in Fig. 9, the forces result from cells that are located in the square-shaped scar region that is enclosed by the red line segments.

Fig. 9

The displacement results are shown, which are solved from and when cells are located according to the cell density function . Hence, implementing algorithm in Fig. 4 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

The pulling forces that are exerted by the cells cause a displacement field over the entire domain of computation (the entire tissue region). The red line segments, which indicate the boundary between scar (containing the cells) and undamaged tissue (not containing cells that exert pulling forces), are displaced as a result of the displacement field caused by the exertion of pulling forces by the cells. The displaced interface between the damaged tissue and scar region is indicated by the black curve enclosed by the red line segments. Since on an average, the cell position is approximately in the center of the domain, the displacements are, on an average, directed towards the center. Furthermore, since the displacement is set to zero at the boundaries of the domain of computation, the magnitude of the displacement vector decreases away from the center. Since the distance between the center and the positions on the red line segments is minimal at the midpoints and maximal at the corners of the red line segments. This causes the contraction and the curved shape of the edge between the scar and undamaged region (the black curve within the red square). In both cases, the finite element solution is in due to the mollified force expressions (gradient of Gaussian distribution), which do not admit any jump discontinuities of the displacement field, and as a result of the use of Lagrangian elements. If the forces exerted by the cells cease, then the boundary between the scar and undamaged region retract back to the red line segments.

Fig. 9 shows the numerical results regarding two approaches. There is no significant difference if the same mesh resolution is used. As the mesh is refined, the solution to the SP approach is smoother, since the “ring” in the center becomes more dominant. In Table 4, it can be concluded that there are no significant differences between the two approaches except for the computational efficiency and the convergence rate of -norm. If the mesh is not fine enough, then the solution to the SP approach is less smooth, hence, the determination of the gradient of the solution is less accurate.

Table 4

Numerical results of two approaches in two dimensions with Gaussian distribution for the positions of biological cells. Figure 3 is implemented and there are 440 biological cells in the computational domain

	SP Approach	Cell Density Approach
L^2-norm of the solution u
‖u^h‖_L^2((0,L))	11.14304909846569	13.188441094735877
‖u^h/2‖_L^2((0,L))	11.180094316889074	13.248171247805358
‖u^h/4‖_L^2((0,L))	11.19130107517498	13.264669660152608
Convergence rate of L^2-norm	1.724918322	1.856132219
‖u‖H^1((0,L))h	12.77795210802095	15.533928099123479
Reduction ratio of the subdomain Ω_w (%)	19.65854	20.55949
Relative ratio of the subdomain Ω_w after deformation (%)	80.34146	79.44051
Time cost (s)	0.62347	0.017315

For Simulation (2), 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 Fig. 4. In Fig. 10, 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 6 displays more details about the two approaches regarding the numerical analysis: most data are more or less the same. 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.

Fig. 10

The displacement results are shown, which are solved from and when cells are randomly located in the subdomain . In other words, it is impossible to write the analytical expression of , subsequently, the algorithm in Fig. 4 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

Table 6

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 Fig. 4 is used

	SP Approach	Cell Density Approach
‖u^h‖_L^2((0,L))	10.105858093727422	12.314518769308366
Reduction ratio of the subdomain Ω_w (%)	24.24192	23.33667
Relative ratio of the subdomain Ω_w (%)	75.75808	76.66333
Time cost (s)	4.20677	1.77219×10^-2

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 and we proved the convergence and the consistency between the solutions to these two approaches in Sect. 2.3.1. In both one and two dimensions, the numerical solutions delivered by finite-element methods with Lagrangian 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.