Mathematical Model and Numerical Simulation for Tissue Growth on Bioscaffolds

Tissue growth on bioscaffolds can be controlled using substrate geometry such as substrate curvature. In this study, we present a mathematical model and numerical simulation method for tissue growth on a bioscaffold to investigate the effect of local curvature on tissue growth. The mathematical model is based on the Allen–Cahn (AC) equation, which has been extensively used to model many problems involving motion by mean curvature. By solving the AC equation using the explicit Euler method, the proposed method is simple and fast. Numerical simulations on various geometries are presented to demonstrate the applicability of the proposed framework on tissue growth on a bioscaffold.


Introduction
In some fracture patients, the bone defects do not get well in the ordinary recuperation period or do not mend at all. In [1], the authors investigated the causes of impaired healing and proposed potential treatment strategies with the combined in silico-in vivo approach. In the field of tissue engineering, bioscaffolds are used as hosts for cell proliferation and extracellular matrix formation. There have been studies on the effects of materials of scaffold substrates [2]. Tissue growth is known to be controlled by substrate geometry such as roughness [3,4] and curvature [5][6][7] of the substrate. Curvature affects cell density because of the contraction or expansion of the local surface area [8].
In previous studies, simulations of tissue growth on scaffold have been performed [9][10][11][12] and applied in many fields [13][14][15][16][17]. Figure 1 shows new tissue formation on various scaffolds (triangular, square, and hexagonal shapes) [18]. Reprinted from Rumpler et al. [18], with permission from The Royal Society. The authors in [18] performed in vitro cell growth studies with osteogenic cells grown on biologically inductive substrate, hydroxylapatite. These cells, pre-osteoblastic mouse MC3T3-E1 cells, exhibit logarithmic growth when cultured in osteogenic conditions. In addition, the authors investigated the geometric features of tissue growth in aspects both the channel surface size and the shape.
The main purpose of this study is to present a mathematical model and numerical simulation for tissue growth on a bioscaffold in two-dimensional (2D) and three-dimensional (3D) space. The mathematical model is based on the Allen-Cahn (AC) equation [19]: where Ω ⊂ R d (d = 2, 3) is a domain with an outer unit normal vector n to ∂Ω. Note that Equation (2) can be replaced with other conditions such as the Dirichlet boundary condition or the periodic boundary condition. Here, φ(x, t) is the difference between the concentrations of the tissue and the ambient medium, µ(x, t) is the chemical potential, F(φ) = 0.25(φ 2 − 1) 2 , and is a positive constant. The AC equation was originally introduced as a mathematical model for antiphase domain coarsening in a binary alloy. The AC equation has been successfully used to model a class of problems such as crystal growth [20,21], image processing [22,23], mean curvature flows [24], and a mixture of two incompressible fluids [25]. We propose a simple mathematical model for tissue growth on scaffolds using features of the AC equation, that is, motion by mean curvature. This paper is organized as follows. In Section 2, we propose the mathematical model for tissue growth on scaffolds. In Section 3, we describe a numerical solution for solving the mathematical model. We present the numerical results for several examples in Section 4. In Section 5, the conclusions are drawn.

Mathematical Model
We briefly review motion by mean curvature of the AC equation [26]. The scalar field φ(x, y, z, t) can be represented by a local coordinate r(x, y, z, t); hence we can represent φ(r) near φ = 0 in terms of the local coordinate r. Figure 2 illustrates the interface between two phases with local coordinate r. Let m = −∇φ/|∇φ| and φ r = ∇φ · m = −|∇φ| be the directional derivative in the direction of the vector m; then, we can rewrite ∆φ as a function of r: where κ = −∇ · m = κ 1 + κ 2 is the mean curvature; and κ 1 and κ 2 are the principal curvatures of the surface. Using this result, Equation (1) becomes When the interfacial transition layer is in the local equilibrium state, then −F (φ)/ 2 + φ rr ≈ 0. Therefore, Equation (3) becomes Let Γ t = {(x, y, z)| φ(x, y, z, t) = 0} be the interface at time t, then its velocity is given by which implies φ r r t + φ t = 0 and therefore r t = −φ t /φ r = κ from Equation (4). From this result, we have V = κ, where V is the normal velocity of the surface. Because φ r ≤ 0 and Equations (1) and (4), both the signs of −µ and κ are the same. We propose the following governing equation for tissue growth under the mean curvature: where µ(x, t) = F (φ(x, t))/ 2 − ∆φ(x, t). That is, if a local mean curvature is positive, then a tissue grows, otherwise it does not [8,18].

Numerical Solution
We present an explicit numerical scheme for the AC equation in the 2D space Ω = (a, b) × (c, d). Performing the numerical computation with the explicit method is efficient because the evolution depends on the a sign of curvature at each grid point. The 3D case is defined analogously. Let N x and N y be positive integers, h = (b − a)/N x = (d − c)/N y be the uniform mesh size, and Because the governing equation is a second-order partial differential equation (PDE), solving the AC equation by the explicit Euler method gives not only accuracy by using a small enough time step size, but also simplicity. For more details about the reason for using the explicit scheme for the second-order PDE, please refer to the article of Jeong and Kim [27]. The advantages of using the scheme are consistent with our purpose for simulation of the tissue growth.
We use two kinds of boundary conditions in this study. One is the Neumann boundary condition, and the other is the periodic boundary condition,

Numerical Experiments
In this section, we perform several numerical tests using various bioscaffold geometries. First, we define the interfacial length parameter m = mh/[2 √ 2 tanh −1 (0.9)] ≈ 0.24015mh, which implies that we have an approximately mh transition layer width [28]. For all tests, we use = m for some integer m, unless otherwise specified. In addition, unless otherwise noted, the zero Neumann boundary condition is used. We note that time-frames for the different simulations will be different because of the increased complexity of the volume component.

Simulation of Tissue Growth on Scaffold In 2D
First, we simulate tissue growth in 2D space to validate the proposed mathematical model although the real scaffolds are 3D. We consider a relatively simple 2D scaffold in the domain Ω = (0, 1) × (0, 1). In all simulations below, the scaffold region is the region in which the value of a scalar field ψ is greater than zero. The simple 2D scaffold is written as follows: which is shown in Figure 3a as the dark region. The surface is initially seeded with tissue layer φ(x, y, 0): which is shown in Figure 3a as the green region. Note that for computational simplicity, we set In Figure 3b,c, we show simulation results of tissue growth on the scaffold in 2D at t = 4000∆t and t = 8000∆t, respectively. Here, the space step h = 1/50, the time step size ∆t = 0.1h 2 , and = 4 are used. The simulation results agree well with the previous solutions [9].
Next, we perform a simulation on a 2D scaffold whose surface is given by a cosine curve which has concavity and convexity to verify the effect of curvature on tissue growth. The 2D scaffold in Ω = (0, 1) × (−0.4, 0.4) is given by The surface is initially seeded with φ(x, y, 0): and the initial cell layer is shown as a green region in Figure 4a. Here, h = 1/200, ∆t = 0.1h 2 , and = 4 are used. As shown in Figure 4, the tissue grows where the curvature is positive. Furthermore, we present a simulation of tissue growth on a star-shaped scaffold which is more complex than the above cases. This simulation is carried out in the domain Ω = (−0.5, 0.5) × (−0.5, 0.5) with h = 1/128, ∆t = 0.1h 2 , and = 4 . An initial star-shaped scaffold is ψ > 0, where ψ(x, y) is given by The initial cell layer is given by where θ is defined as above. Figure 5 shows the temporal evolution of tissue on the star-shaped scaffold.
(a) (b) (c) Figure 5. (a-c) are snapshots of tissue growth simulation with the star-shaped scaffold in 2D at t = 0, 600∆t, and 1200∆t, respectively. The green region is the tissue layer on the given scaffold which is represented by the dark region. For interpretation of the references to color in the figure, the reader is referred to the web version of the article.

Simulation of Tissue Growth on Scaffold In 3D
We consider tissue growth on a 3D scaffold in Ω = (0, 1) 3 . The numerical scheme is similar to the 2D case. The parameters used are h = 1/50, ∆t = 0.1h 2 , and = 4 . The solid scaffold structure is given as follows: and is shown as a volume enclosed by the dark color surface in Figure 6a. The initial tissue layer is given by and is shown as a volume enclosed by the green colored surface on the scaffold in Figure 6a. Snapshots of tissue growth simulation are shown in Figure 6b,c at t = 3200∆t and 6400∆t, respectively. The tissue growth on the scaffold according to the curvature shows the similar simulation results from the previous test [9].
The simulation is performed in Ω = (−0.5, 0.5) 3 with h = 1/64, ∆t = 0.1h 2 , and = 4 . Furthermore, N x = N y = N z = 64 are used. We use the periodic boundary condition. Figure 7a-c show snapshots of the tissue growth simulation in 3D on the P, D, and G surface scaffolds (from top to bottom) at t = 0, 500∆t, and 1000∆t, respectively.
The temporal evolutions of the volume and those of the surface area, V and S, of cell layer over TPMSs are shown in Figure 8a,b, respectively. The volume and the surface area are defined as respectively. Note that 0.5 is the volume of each scaffold. Here, a smoothed Dirac delta function [30], where |∇ h φ| is defined as follows: and φ y , φ z are similarly defined.
(a) (b) (c) Figure 7. (a-c) are snapshots of tissue growth simulation in 3D on the P, D, and G surface scaffolds (from top to bottom) at t = 0, t = 500∆t, and t = 1000∆t, respectively. The green region is the tissue layer on the given scaffold which is represented by the dark region. For interpretation of the references to color in the figure, the reader is referred to the web version of the article.
Furthermore, the temporal evolution of the averaged local curvature, κ avg , is shown in Figure 9. The averaged local curvature is defined as From the results in Figures 8 and 9, we can observe that scaffolds with initially larger surface areas grow faster than scaffolds with smaller surface areas. Moreover, scaffolds with large surface areas have the larger value of the averaged local curvature than that of scaffolds with small surface areas. We consider the CPU time in seconds for tissue growth on TPMSs. The CPU time taken for the domain to fill up takes 131 s on P surface, 47 s on D surface, and 72 s on G surface. This CPU time results are consistent with the results in Figure 8a.

Comparison between Experiments and Simulations
In this section, we conduct a comparative study with the actual experimental data to validate our numerical simulation. There is an actual experimental study on whether geometric features affect tissue growth by using four different cross-sectional shapes (triangle, square, hexagon, and circle) in [18]. Figure 10 illustrates the domain of experiment which consists of initial perimeter P and projected tissue area A. Because we do consider the effects of chemicals, we omit the description of substrate. The experiment is conducted by the assumption that any position of the tissue surface normally moves to the surface by distance ds during the time step dt. A simple 2D mathematical formulation is expressed as follows in [18]: ds/dt = −λκ, where λ is a constant of proportionality and κ is local curvature. We conduct a numerical simulation under the same conditions to match ours with the experiment, i.e., λ = 0.97 and P = 3.14 are taken for the numerical test. Note that we conduct the simulation only for P = 3.14 for simplicity of exposition. As shown in Figure 11, the average thickness of cell, A/P, is good agreement with the experimental result. Figure 10. P is the initial perimeter of the channel cross-section and A is the projected tissue area. Reprinted from Rumpler et al. [18], with permission from The Royal Society. Figure 11. Comparison of numerical simulation and experiment in [18] with P = 3.14. Four different channel shapes (square, hexagon, triangle, and circle) are represented by colored marker-line (magenta, blue, green, and red) from top to bottom, respectively. Reprinted from Rumpler et al. [18] with permission from The Royal Society.

Conclusions
We proposed a mathematical model for tissue growth on a bioscaffold in 2D and 3D. The AC equation which has motion by mean curvature is a major base for this modeling. Because we solve the AC equation by the explicit Euler method, the proposed method is simple and accurate. We showed several numerical simulations in 2D and 3D by using our proposed modeling. Tissue on scaffolds grew according to the sign of curvature, and these simulations agree well with previous biological studies. In particular, simulations of tissue growth on TPMSs provide interesting results, and we investigate which kind of TPMSs is the most suitable for tissue growth on a bioscaffold. Therefore, we present a mathematical framework of tissue growth on various bioscaffolds using the AC equation. In future work, we will aim to model the chemistry on the substrate because the chemistry of the surface is important for the growth of a tissue on scaffold.
Author Contributions: All authors, H.G.L., J.P., S.Y., C.L., and J.K., contributed equally to this work and critically reviewed the manuscript. Acknowledgments: The authors thank the reviewers for their constructive and helpful comments on the revision of this article.

Conflicts of Interest:
The authors declare no conflict of interest.