6 Computational modeling of the neotissue growth process‣ Mathematical and computational modeling of metallic biomaterials biodegradation

6.1 Introduction

Neotissue is defined as the cells and the extracellular matrix they produce. Coupling biodegradation and neotissue growth models can be useful for tuning the biodegradation rate to the rate of regeneration of new tissue, an example of which can be found in Chapter 10. Similar to the biodegradation model, the modeling of the tissue formation process can take advantage of the interface tracking techniques in which the neotissue surface evolves over time, representing the growth process. Neotissue growth in porous scaffolds has been shown to be depending on the local mean curvature of the interface between the neotissue and the surrounding void space [36, 37, 221]. The commonly-used interface tracking techniques provide an efficient way to formulate curvature-driven problems. Although these models do not have immediate clinical applications, they can be quite useful for in vitro tissue engineering experiments where the proliferation behavior and growth of various cell types need to be improved.

Numerical tracking of interface movement has been widely used for certain modeling applications in science and engineering for multi-material and multiphase problems such as solidification, melting, corrosion, and grain growth to name a few [251]. The most popular Eulerian methodologies in this regard are the level-set [203, 16, 219], volume-of-fluid (VOF) [218], and phase-field methods [40, 33].

The basic idea of the level-set method is to employ the Hamilton-Jacobi (HJ) algorithm for solving the general interface advection equation. The independent variable in the level-set method is a signed distance function called the level-set function $\psi$ [219]. The level-set function should be re-initialized as the interface evolves, which is the reason behind an inadvertent mass loss, one of the most prominent shortcomings of the level-set method. Although the VOF method is not vulnerable to the mass loss issue, calculating the interface curvature is difficult from the volume fraction [251], making it less efficient to be used for curvature-driven problems.

To overcome these challenges, diffuse interface methods [15] have gained attention in recent decades, among which the phase-field method has shown potential for solving complex moving interface problems. Contrary to the level-set method, the interface is considered a smooth transition between phases, which usually has a finite width in diffuse interface methods. In the phase-field method, a non-conserved (or conserved) order parameter $\phi$ is defined such that $\phi=1$ in one bulk phase and $\phi=-1$ in the other. Then, the smooth transition between these two phases ( $-1<\phi<1$ ) is marked as the interface. One of the advantages of the phase-field method is that the derived equation can be solved over the entire desired domain without considering the location of the interface. Moreover, although the curvature and interface normal vectors are not formulated explicitly, the phase-field method is suitable for problems in which the evolution of the interface depends on the local curvature or a field acting normal to the interface [251]. The phase-field method keeps a constant thickness for the smooth transition region normal to the interface, and as a result, no re-initialization as for the level-set method is needed.

The phase-field method has been already proved to be an efficient interface tracking technique for various problems in micro/meso scales such as solidification [143, 40], microstructural evolution [51], grain growth [50], crack propagation [114, 246], electromigration [35] and extractive metallurgy [33]. However, it has been used recently for dealing with problems described in macro level, such as corrosion [183, 168, 126, 169, 17, 255, 47] and cell/tissue growth [135, 158]. This chapter discusses the development of a phase-field model of the tissue growth process to describe the cell growth behavior on 3D surfaces as a moving-boundary problem. Additionally, a similar model was developed based on the level-set method to compare the performance and results of both interface tracking methods. Both models were implemented using the finite element method.

6.2 Deriving the model

This section demonstrates the derivation of the phase-field and level-set equations from the general advection equation, showing the similarities and differences of these interface tracking techniques for moving boundary problems.

6.2.1 General equation of interface motion

The general interface advection equation for an Eulerian description of interface movement can be written as [251]:

\frac{\partial\phi}{\partial t}+\boldsymbol{u}\cdot\nabla\phi=0,

(6.1)

where $\phi$ is the phase-field and $\boldsymbol{u}$ is the interface velocity. The velocity $\boldsymbol{u}$ can be divided into normal ( $u_{\mathrm{n}}$ ) and external velocity components ( $\boldsymbol{u}_{\mathrm{e}}$ ):

\boldsymbol{u}=u_{\mathrm{n}}\boldsymbol{n}+\boldsymbol{u}_{\mathrm{e}},

(6.2)

in which $\boldsymbol{n}=\nabla\phi/|\nabla\phi|$ is the unit vector normal to the interface. So, Eq. 6.1 can be rewritten as:

\frac{\partial\phi}{\partial t}+u_{\mathrm{n}}|\nabla\phi|+\boldsymbol{u}_{% \mathrm{e}}\cdot\nabla\phi=0.

(6.3)

The normal velocity can be decomposed into more components to take into account the effect of interface curvature ( $\kappa$ ) such that the terms are independent and proportional to the curvature, respectively:

u_{\mathrm{n}}=a-b\kappa,

(6.4)

where the coefficients $a$ and $b$ have units of $\mathrm{m}/\mathrm{s}$ and $\mathrm{m}^{2}/\mathrm{s}$ . Substituting this into Eq. 6.2 yields the final form of the interface motion equation:

\frac{\partial\phi}{\partial t}+a|\nabla\phi|+\boldsymbol{u}_{\mathrm{e}}\cdot% \nabla\phi=b\kappa|\nabla\phi|.

(6.5)

6.2.2 Phase-field formulation

To further proceed with the phase-field formulation, a proper kernel should be selected for the phase-field variable, which can be done based on Beckermann et al. [32]:

\phi=-\tanh\left(\frac{n}{\sqrt{2}w}\right),

(6.6)

in which $w$ is the thickness of the transition profile ( $\phi$ varies from $-0.9$ to $+0.9$ in a narrow layer with the width of $3\sqrt{2}w$ ), and $n$ is the coordinate normal to the interface. The curvature can be written as a function of the phase-field variable:

\kappa=\nabla\cdot\boldsymbol{n}=\nabla\cdot\left(\frac{\nabla\phi}{|\nabla% \phi|}\right)=\frac{1}{|\nabla\phi|}\left[\nabla^{2}\phi-\frac{(\nabla\phi% \cdot\nabla)|\nabla\phi|}{|\nabla\phi|}\right].

(6.7)

Using the defined kernel, the terms in Eq. 6.7 can be expressed as:

|\nabla\phi|=-\frac{\partial\phi}{\partial n}=\frac{1-\phi^{2}}{\sqrt{2}w}% \quad\text{ and }\quad\frac{(\nabla\phi\cdot\nabla)|\nabla\phi|}{|\nabla\phi|}% =\frac{\partial^{2}\phi}{\partial n^{2}}=-\frac{\phi\left(1-\phi^{2}\right)}{w% ^{2}}.

(6.8)

Substituting Eq. 6.8 into Eq. 6.7 yields to the following definition of interface curvature:

\kappa=\frac{1}{|\nabla\phi|}\left[\nabla^{2}\phi+\frac{\phi\left(1-\phi^{2}% \right)}{w^{2}}\right],

(6.9)

which subsequently changes Eq. 6.5 into:

\frac{\partial\phi}{\partial t}+a|\nabla\phi|+\boldsymbol{u}_{\mathrm{e}}\cdot% \nabla\phi=b\left[\nabla^{2}\phi+\frac{\phi\left(1-\phi^{2}\right)}{w^{2}}% \right].

(6.10)

Eq. 6.10 is the derived form of the phase-field equation for tracking an evolving interface, containing terms corresponding to normal interface motion, advection by an external field, and curvature-driven movement. The term $|\nabla\phi|$ in Eq. 6.10 can be replaced by its definition in Eq. 6.8 to form another version of the equation:

\frac{\partial\phi}{\partial t}+a\frac{1-\phi^{2}}{\sqrt{2}w}+\boldsymbol{u}_{% \mathrm{e}}\cdot\nabla\phi=b\left[\nabla^{2}\phi+\frac{\phi\left(1-\phi^{2}% \right)}{w^{2}}\right],

(6.11)

which is an easier version to be implemented using numerical methods. The unique term on the right-hand side of Eq. 6.11 is a characteristic of the phase-field method.

6.2.3 Level-set formulation

From the mathematical perspective, the level-set equation has a direct connection to the phase-field equation and can be derived by replacing the phase-field variable with a sign distance function. To this end, Eq. 6.5 can be rewritten to be a level-set equation:

\frac{\partial\psi}{\partial t}+a|\nabla\psi|+\boldsymbol{u}_{\mathrm{e}}\cdot% \nabla\psi=b\kappa|\nabla\psi|,

(6.12)

with $\psi$ being a sign distance function that descibes the distance of each point of the computational domain to the interface. This implies that the zero iso-contour of the level-set function defines the interface.

6.3 Dimensionless forms for various cases

6.3.1 Stationary interface

A stationary interface, where there is no interface motion ( $a=0$ and $\boldsymbol{u}_{\mathrm{e}}=0$ ), is a primary problem in examining the formulation and select proper grid spacing parameters. For a 1-D case, Eq. 6.11 can be simplified as:

\frac{\partial\phi}{\partial t}=b\left(\frac{\partial^{2}\phi}{\partial x^{2}}% +\frac{\phi\left(1-\phi^{2}\right)}{w^{2}}\right).

(6.13)

To scale this equation, the following dimensionless variables can be defined:

x^{\prime}=\frac{x}{x_{c}}\quad\text{and}\quad t^{\prime}=\frac{t}{t_{c}}.

(6.14)

So, Eq. 6.13 can be rewritten using these new variables:

\frac{1}{t_{c}}\frac{\partial\phi}{\partial t^{\prime}}=\frac{b}{x_{c}^{2}}% \frac{\partial^{2}\phi}{\partial x^{\prime 2}}+\frac{b}{w^{2}}f(\phi),

(6.15)

with $f(\phi)$ being defined as:

f(\phi)=\phi\left(1-\phi^{2}\right).

(6.16)

Defining $x_{c}=w$ and $t_{c}=w^{2}/b$ leads to the following dimensionless form of Eq. 6.13:

\frac{\partial\phi}{\partial t^{\prime}}=\frac{\partial^{2}\phi}{\partial x^{% \prime 2}}+\phi\left(1-\phi^{2}\right).

(6.17)

Numerical results for Eq. 6.17 is depicted in Fig. 6.1, and the phase-field profile is compared with the level-set distance function profile. An appropriate value for grid spacing and the layer width should be selected such that $0.25w<\Delta x^{\prime}<0.5w$ . Additionally, the selected value of $w$ should satisfy $w<R/4.2$ , in which $R$ is the local radius of curvature [251].

Figure 6.1: Comparison of phase-field variable and level-set function on the interface of a stationary interface with step function like initial condition [251]

6.3.2 Evolution under constant normal speed

For a problem in which the interface moves with constant velocity exclusively, Eq. 6.11 in 1-D can be simplified according to the condition of $a=\text{const.}$ , $\boldsymbol{u}_{\mathrm{e}}=0$ , and $b=0$ :

\frac{\partial\phi}{\partial t}+a\frac{1-\phi^{2}}{\sqrt{2}w}=\beta\frac{% \partial^{2}\phi}{\partial x}+\beta\frac{f(\phi)}{w^{2}},

(6.18)

in which $\beta$ is a numerical parameter for smoothing the interface and relaxation behavior of the phase-field profile.

Defining dimensionless variables according to Eq. 6.14 yields to:

\frac{1}{t_{c}}\frac{\partial\phi}{\partial t^{\prime}}+\frac{a}{w}\frac{1-% \phi^{2}}{\sqrt{2}}=\frac{\beta}{w^{2}}\frac{\partial^{2}\phi}{\partial x^{% \prime 2}}+\frac{\beta}{w^{2}}f(\phi),

(6.19)

which can be reordered to:

\frac{\partial\phi}{\partial t^{\prime}}+\frac{1-\phi^{2}}{\sqrt{2}}=\beta^{% \prime}\left(\frac{\partial^{2}\phi}{\partial x^{\prime 2}}+f(\phi)\right),

(6.20)

with $\beta^{\prime}=\beta/aw$ . For a stable numerical implementation, $\Delta t^{\prime}/\Delta x^{\prime}<0.1$ and $\beta^{\prime}<1.2$ should met roughly [251].

6.3.3 Curvature-driven interface evolution

A curvature-driven motion, which is desired for the current study, is straightforward to formulate using the phase-field method. A dimensionless form of Eq. 6.11 can be derived using a similar method for the stationary interface for a multidimensional case with $u_{\mathrm{n}}=-b\kappa$ , $a=0$ , and $\boldsymbol{u}_{\mathrm{e}}=0$ :

\frac{\partial\phi}{\partial t^{\prime}}=\nabla^{\prime 2}\phi+\phi\left(1-% \phi^{2}\right),

(6.21)

in which $t^{\prime}$ and $\nabla^{\prime}$ are defined similar to Eq. 6.14 as $t^{\prime}=t/(w^{2}/b)$ and $\nabla^{\prime}=\nabla/w$ , respectively.

6.4 Adapting the formulation for curvature-driven tissue growth

Due to intrinsic support of interface curvature in the phase-field and level-set methods, an in silico model of curvature-based tissue growth can be efficiently implemented using these principles. The growth-induced changes in the neotissue topology during the culture process can be seen as a moving interface between two different domains [221]. In this study, one domain represents the neotissue volume ( $\Omega_{\text{nt}}$ ), and the other one is the void ( $\Omega_{\text{v}}$ ), which are separated by an interface ( $\Gamma$ ) as can be seen in Fig. 6.2.

Figure 6.2: Schematic representation of the phase-field and level-set models for tissue growth, in which the neotissue domain ( $\Omega_{\text{nt}}$ ) is separated from the void domain ( $\Omega_{\text{v}}$ ) by an interface ( $\Gamma$ ).

The interface $\Gamma$ evolves over time to fill the void space, having a faster growth in regions with higher curvature. Based on this definition, the phase-field variable can be defined as follows to separate these domains:

\left\{\begin{array}[]{ll}\phi=1&\text{ in }\Omega_{\text{nt}}\\ \phi=-1&\text{ in }\Omega_{\text{v}}\\ -1<\phi<1&\text{ in }\Gamma\end{array}\right.

(6.22)

Similarly, a level-set function can be defined such that it separates the neotissue and void domains:

\left\{\begin{array}[]{ll}\psi>0&\text{ in }\Omega_{\text{nt}}\\ \psi<0&\text{ in }\Omega_{\text{v}}\\ \psi=0&\text{ in }\Gamma\end{array}\right.

(6.23)

In order to adapt Eq. 6.21 for the curvature-driven process of neotissue growth, one can consider the following equation:

\frac{\partial\phi}{\partial t^{\prime}}=\left(\nabla^{\prime 2}\phi+\phi\left% (1-\phi^{2}\right)\right).H\left(\nabla^{\prime 2}\phi+\phi\left(1-\phi^{2}% \right)>0\right),

(6.24)

in which $H$ denotes a Heaviside step function. Eq. 6.24 implies that the growth is only allowed for regions with a positive curvature (right hand side of Eqs. 6.11 and 6.21).

Using the same approach, a similar level-set formulation can be obtained based on Eq. 6.12 by omitting the normal velocity and curvature terms and embedding the effect of the curvature in the velocity field. Doing this yields a convection equation for the distance function:

\frac{\partial\psi}{\partial t}+\boldsymbol{u}\cdot\nabla\psi=0,

(6.25)

in which the convection velocity field can be defined as:

\boldsymbol{u}=\left\{\begin{array}[]{ll}-\kappa\boldsymbol{n}&\text{ if }% \kappa>0\\ 0&\text{ if }\kappa\leq 0\end{array}\right.

(6.26)

with $\kappa$ being calculated similarly to Eq. 6.7 for a distance function $\psi$ . This implies that neotissue grows faster where the curvature is higher and does not grow if the curvature is negative or equal to zero [36]. The negative sign in Eq. 6.26 is due to our definition of $\psi$ , where the normal $n_{\Gamma}$ points toward neotissue, so growth has to be towards the opposite of the gradient of the level-set function ( $\nabla\psi$ ).

6.5 Numerical implementation

6.5.1 Phase-field model

The numerical solution of the phase-field equation requires dealing with the nonlinearity of the equation. Additionally, in the case of the dimensional form (Eq. 6.11), small coefficients of the state variable in the PDE lead to numerical difficulties. As a result, numerical implementation of the phase-field equation, especially for the spectral methods and the finite element method, is not straightforward and is an active field of research [235, 2].

In the finite element method, the solution of a PDE is calculated based on a sum of a set of certain basis functions, which are commonly piecewise linear, quadratic or polynomial functions that are non-zero only on a small element. For doing this, the PDE is first written in a weak formulation, and then the weak form is projected on a discretized space (a set of elements) to be written as the summation of the basis functions.

In this section, the numerical solution of the stationary form (Eq. 6.13) and its corresponding considerations are elaborated as an example of employing the finite element formulation for simulating the phase-field equation. So, by assuming $b=1$ , the problem can be summarized as:

\left\{\begin{array}[]{ll}\frac{\partial\phi}{\partial t}-\Delta\phi+\frac{1}{% w^{2}}f(\phi)=0,&(x,t)\in\Omega\times(0,T]\\ \left.\frac{\partial\phi}{\partial n}\right|_{\partial\Omega}=0\\ \left.\phi\right|_{t=0}=\phi_{0}\end{array}\right.

(6.27)

which demonstrates the PDE, the boundary condition, and the initial condition of the phase-field variable where $\Omega$ is the domain of interest, $\partial\Omega$ is its boundary, and $T$ is the final time. Deriving the weak formulation of Eq. 6.27 is relatively straightforward as it can be seen as a time-dependent diffusion-reaction PDE, but the difficulty arises in choosing the numerical stability scheme for discretizing the temporal derivative and dealing with the nonlinearity of $f(\phi)$ when normally the $\frac{1}{w^{2}}$ coefficient is a small number.

Incorporating a first-order semi-explicit scheme for Eq. 6.27 yields [2]:

\frac{1}{\Delta t}\left(\phi^{n+1}-\phi^{n},v\right)+\left(\nabla\phi^{n+1},% \nabla v\right)+\frac{1}{w^{2}}\left(f\left(\phi^{n}\right),v\right)=0,\quad% \forall v\in H^{1}(\Omega),

(6.28)

where $\Delta t$ is the time step, $(\cdot,\cdot)$ denotes the inner product, and $H^{1}(\Omega)$ is the Sobolev space of the domain $\Omega$ , which is a space of functions whose derivatives are square-integrable functions in $\Omega$ . The main issue with this discretization scheme is its restrictive time step condition which should satisfy [235]:

\Delta t<\frac{2w^{2}}{L},

(6.29)

where $L$ is a limit related to the non-linear part:

\max\left|f^{\prime}(\phi)\right|\leq L.

(6.30)

Obviously, since $\Delta t\sim w^{2}$ , a very small time step is required to achieve stability in this scheme.

Taking advantage of a fully implicit scheme improves the stability because it will be unconditionally stable, but it results in an equation that is difficult to implement as it needs to solve a fixed point problem at each time step. For example, a modified second-order implicit Crank-Nicolson scheme for Eq. 6.27 can be written as [2, 73]:

\left(\frac{\phi^{n+1}-\phi^{n}}{\Delta t},v\right)+\left(\nabla\frac{\phi^{n+% 1}+\phi^{n}}{2},\nabla v\right)+\frac{1}{w^{2}}\left(\tilde{f}\left(\phi^{n+1}% ,\phi^{n}\right),v\right)=0,\quad\forall v\in H^{1},

(6.31)

where:

\tilde{f}(u,v)=\left\{\begin{array}[]{ll}\frac{F(u)-F(v)}{u-v}&\text{ if }u% \neq v\\ f(u)&\text{ if }u=v\end{array}\right.

(6.32)

in which $F$ is the potential term ( $f(\phi)=F^{\prime}(\phi)$ ).

An alternative can be deriving a stabilized semi-implicit scheme by adding a stabilization term to Eq. 6.28. The first-order version of such a scheme can be written as:

\left(\frac{1}{\Delta t}+\frac{S}{w^{2}}\right)\left(\phi^{n+1}-\phi^{n},v% \right)+\left(\nabla\phi^{n+1},\nabla v\right)+\frac{1}{w^{2}}\left(f\left(% \phi^{n}\right),v\right)=0,\quad\forall v\in H^{1}(\Omega),

(6.33)

which is unconditionally stable for any $S\geq\frac{L}{2}$ [235].

6.5.2 Level-set model

The derived level-set PDE (Eq. 6.25) is an advection equation, which can be implemented numerically using the finite element method, in which the temporal term is discretized by the backward Euler method, and the advection term can be treated with the method of characteristics.

A key parameter of the developed model is the local growth velocity of the neotissue. In the current implementation, the velocity was dependent on the interface’s local mean curvature as shown in [37, 101]. In order to match the growth velocity to experimental data, a coefficient can be added to the derived interface convection velocity (Eq. 6.26):

\boldsymbol{u}=\left\{\begin{array}[]{ll}-\kappa A\boldsymbol{n}&\text{ if }% \kappa>0\\ 0&\text{ if }\kappa\leq 0\end{array}\right.

(6.34)

The model calibration performed by Guyot et al. [101] was for a culture condition on titanium scaffolds in a bioreactor setting, estimating parameter $A$ to be $4\times 10^{-14}\text{m}^{2}/\text{s}$ , obtained using trial and error from the experimental data on low flow rate tests [205]. More dedicated calibration experiments were performed on prismatic structures, demonstrating a considerably slower growth on the CaP scaffolds, nevertheless confirming the curvature-based nature of tissue growth [113].

In practical implementations, the distance function is not differentiable at every location of the domain due to discontinuities in the gradients, so one can consider taking advantage of artificial diffusion terms to overcome this issue, leading to the following equations for the normal vector and curvature calculation:

\boldsymbol{n}=\frac{\nabla\varphi}{|\nabla\varphi|}+\varepsilon\Delta% \boldsymbol{n}

(6.35)

\kappa=\nabla\cdot\boldsymbol{n}+\varepsilon\Delta\kappa,

(6.36)

in which $\varepsilon$ denotes the numerical diffusion coefficient.

6.6 Simulation setup

Numerical simulations of neotissue formation on various shapes (scaffolds) were carried out using the developed phase-field and level-set models to compare their performance. In 2D, two shapes were used, a square and a semi-circle (Fig. 6.3), which were simulated using both phase-field and level-set models. The qualitative comparison of the two developed models were done using this 2D setup, so for 3D cases, the models were simulated on vastly different geometries for just checking the performance of the interface tracking techniques on 3D scaffolds. In 3D, both models were evaluated on a cube, but the level-set model was also used to simulate the cell growth behavior on scaffolds with triply periodic minimal surface (TPMS) shapes.

Figure 6.3: Schematic representation of the simulation domains

The initial configuration of the phase-field variable and level-set distance function corresponds to a homogenous single cell layer over the scaffold struts with a thickness equal to $10\mu m$ in a dimensional setup [64]. This is depicted as the green layer in Fig. 6.3. For the 3D level-set model, neotissue growth was simulated for a variety of TPMS family and compared in a qualitative manner. A full quantitative prediction is not possible due to the absence of relevant validation experiments, which explains why comparisons between geometries are made over non-dimensional time.

The derived weak forms were implemented using FreeFEM open-source PDE solver [112]. An Eulerian computational mesh for each simulation was constructed by generating tetrahedral elements using the internal mesh generator of FreeFEM, called BAMG. To decrease run time and increase the performance of the simulation, the mesh was partitioned using METIS graph partitioner [144] and HPDDM preconditioner [140], available in FreeFEM and the PETSc toolkit [25]. Moreover, the efficiency was boosted by using the HYPRE BoomerAMG preconditioner [78] and the GMRES iterative solver [223] of the PETSc toolkit.

6.7 Results and discussion

In the current study, various 2D and 3D cases for the simulation of curvature-based neotissue formation were prepared and separately simulated using the phase-field and level-set models. The cases were a square and a semi-circle for 2D simulations, mimicking the situation in which cells were homogeneously seeded onto the scaffold to initiate the production of extracellular matrix. Similarly, a 3D cube was considered for evaluating the performance of the implementation of both models in 3D. Additionally, the level-set model was used to simulate cell proliferation on TMPS scaffolds with quantitative measurements being carried out for the percentage and filling rate of generation of neotissue.

The results of neotissue growth model simulations can be visualized by depicting the +1 part of the phase-field variable (Eq. 6.22) and the positive part of the level-set function (Eq. 6.23). Fig. 6.4 demonstrates such a visualization for the phase-field and level-set simulations performed on the square domain. Cells were seeded on the perimeter of the square, and the formation of neotissue was modeled using the evolution of the phase-field variable and level-set function. The light gray region in these figures shows the +1 part of the phase-field variable and the positive part of the level-set function.

Figure 6.4: Simulation result of the phase-field (A) and level-set (B) models for neotissue growth in the square domain. The light gray region shows the +1 part of the phase-field variable and the positive part of the distance function for the phase-field (A) and level-set (B) models, respectively. The evolution occurs from left to right over non-dimensionalized time.

As shown in Fig. 6.4, a qualitative comparison of the evolution of the formed tissue interface between the phase-field and level-set models indicates that they capture the curvature-driven growth similarly. It seems that the level-set model resulted in sharper interfaces, in which the surfaces without curvature do not move at all. But, in the phase-field predictions, the surfaces not having any curvature slightly move. This can be seen by comparing the growth pattern on the middle part of the top edge of the square, where it starts to grow only in the fourth column in Fig. 6.4 for the level-set results, while a minor move can be seen in the phase-field predictions. This slightly different behavior can be related to the implementation details, such as the lack of a proper Heaviside function in FreeFEM that the phase-field model depends on. Another possible reason for this difference can be the visualization aspect since the transition region of the phase-field variable representing the surface is thicker in the phase-field model in comparison to the zero iso-contour of the level-set function. As a result, the evolution can be artificially plotted in the visualization software, which was ParaView in this case. This effect can be investigated further by refining the mesh to minimize artificial evolution.

A similar visualization is depicted in Fig. 6.5 for the neotissue formation on the semi-circle domain. Results show that the phase-field and level-set models have good agreement on the way they treat the curvature-based tissue growth. As seen in the figure, both models show identical predictions on the surface with curvature (the curvy edge of the semi-circle). However, the small difference observed for the square case can also be seen here on the top surface, where the neotissue grows slightly faster in the phase-field model.

Figure 6.5: Simulation result of the phase-field (A) and the level-set (B) models for neotissue growth in the semi-circle domain. The evolution occurs from left to right over non-dimensionalized time.

Fig. 6.6 shows the evolution of the phase-field variable corresponding to the results depicted in Fig. 6.5-A, demonstrating how the diffuse interface model works. In this figure, the concept of the diffusive interface can be observed, where a narrow region between the two phases (the neotissue and the void space in this case) is formed and moves over time. This visualization demonstrates the internal mechanism of the phase-field model, in which the narrow region is kept at a fixed length (defined by the kernel in Eq. 6.6 for a length of $3\sqrt{2}w$ ) and gets advected over time by the phase-field equation (Eq. 6.11). In contrast, in the level-set formalism, the interface is not tracked as a fixed length region, and instead, the zero iso-contour of a signed distance function is the interface between the phases. The signed distance function maps each node of the whole space (the desired computational domain) into the distance to the interface, and as a result, the change in the function determines the movement of the interface.

Figure 6.6: Visualization of the evolution of phase-field variable in the semi-circle domain. The evolution is depicted over time in an arbitrary unit, occurring from top to bottom and from left to right in each row. The red, blue, and white show the +1 region, -1 region, and the transient phase representing the interface, respectively.

Similar results are obtained in 3D, where the cells are seeded on the circumstance of a cuboid. Fig. 6.7 shows the simulation results of this 3D case, in which the evolution of the phase-field variable is converted to a bulk of the formed neotissue. As can be seen in this figure, the phase-field model performs well in 3D, showing a similar growth behavior to the 2D model. The mesh in the 3D case is relatively coarse, but the captured behavior of tissue growth by the phase-field model seems to be more acceptable than in the 2D case. This claim can be observed on the top and bottom surfaces of the cube, where the interface does not move until a curvature is created in those regions, a behavior that was captured better by the level-set model in the 2D cases.

Figure 6.7: Simulation result of the phase-field model for neotissue growth in a 3D cube domain, in which the +1 regions of the phase-field variable are visualized to show the formation of neotissue over non-dimensionalized time.

Similarly, Fig. 6.8 shows the evolution of neotissue formation in a 3D level-set simulation, in which the results are similar to the ones obtained using the phase-field model (Fig. 6.7). Generally speaking, the level-set model requires finer mesh in 3D in order to capture the interface movement accurately, which is the reason behind having a smoother interface in Fig. 6.8 in comparison to Fig. 6.7. However, a qualitative comparison shows that both models act equivalently for capturing the curvature-dependent growth in 3D.

Figure 6.8: Simulation result of the level-set model for neotissue growth in a 3D cube domain, in which the positive regions of the level-set variable are visualized to show the formation of neotissue over non-dimensionalized time.

Fig. 6.9 depicts various simulation results of the netissue formation in gyroid-TPMS scaffolds, demonstrating an example of the final application of these models in action. The reason for choosing gyroid-TPMS scaffolds is their promising performance for neotissue formation due to their favourable local curvature. Various geometries for investigating the pore size and wall thickness were evaluated in order to obtain the combination leading to optimal neotissue growth. Pore size was varied between $700\mu m$ and $1.3mm$ , and wall thickness was varied from $200\mu m$ to $800\mu m$ , which are ranges that take into account restrictions of the manufacturing process in terms of the smallest feature dimensions. Balancing the need for swift neotissue ingrowth (Fig. 6.9-B, relevant for short-term implant stability) with the volume of neotissue formed (Fig. 6.9-A, relevant for long-term dental implant stability), the combination of $700\mu m$ pore size and $200\mu m$ wall thickness (Fig. 6.9-C) seemed to be the most optimal structure.

Figure 6.9: Neotissue formation quantified in absolute volume (A) and filling percentage (B) for cylindrical test samples (diameter 6 mm, height 6 mm). The different combinations are indicated in the legend by 2 numbers, the first of which refers to pore size (7: 700 $\mu$ m, 10: 1 mm, 13: 1.3 mm), and the second refers to wall thickness (2: 200 $\mu$ m; 5: 500 $\mu$ m; 8: 800 $\mu$ m;). (C) Side view and cross-sectional view of neotissue growth in gyroid (7-2) scaffold for different levels of filling, starting with the initial condition at t=0 (top). Scale bars: 6 mm.

The study shown as an example above [113] demonstrates the relevance of the developed tissue growth models in tissue engineering applications. In silico modeling is widely used for tissue engineering as it offers a more exhaustive approach compared to a “trial-and-error” method and reduces the number of experimental tests. Optimization of scaffold structures for bone tissue engineering purposes is often corroborated by comparison with in vitro tests [11, 220] and only a small number of in vivo studies have been reported in this regard [173, 175, 176]. However, in these in vivo studies, optimization was first performed on mechanical properties rather than the structural elements such as local curvature underlying the in silico model presented in this chapter. Having such a model developed makes it easier to investigate the structural elements of scaffolds for improved regenerative performance. Using the modeling results makes it possible to limit the number of conditions tested in vivo.

6.8 Challenges in coupling tissue growth and biodegradation models

Coupling our biodegradation model [31] (Chapter 3) with the developed tissue growth model in this chapter could have been a milestone of the current PhD thesis, where the neotissue is produced by the cells seeded on a biodegradable scaffold. Another interesting example could be using biodegradable porous implants in bone healing applications (similar to the model presented in Chapter 10), in which the implant degrades while new bone forms and replaces the porous implant. However, despite the effort put on tissue growth modeling, such coupling can be complicated to accomplish in a mechanistic manner due to specific technical difficulties. In this section, these challenges are briefly reviewed.

Both the biodegradation and the tissue growth models are free boundary problems, in which an interface capturing method (phase-field or level-set for the tissue growth and level-set for the biodegradation) was used to track the movement of a boundary. The moving boundary in the biodegradation model is the corrosion front, the interface between the metallic part and the surrounding environment, which can be a static electrolyte in immersion tests or circulated solutions in a perfusion setup. Various mass transfer boundary conditions are defined on this boundary [31], meaning that the boundary conditions are not fixed and move with the evolution of the surface. In other words, the boundary conditions of a set of equations are defined on the solution of one of the governing equations of the model (level-set equation), which is one of the challenging parts of the implementation of that model. On the other hand, the moving boundary in the tissue growth model is the surface of the formed neotissue, the interface between the neotissue and the surrounding environment in the bioreactor, which is considered as the void space since there is no tissue in that region. Coupling these two models requires defining the moving boundary problem of neotissue growth on the solution of the interface capturing of the biodegradation model. Regardless of using the same technique for both problems or using separate methods for each problem, such coupling can be quite challenging from the implementation perspective.

The second problem is related to the behavior of neotissue formation after coupling the models when the biodegradable material shrinks and a new void space appears at the back of the formed neotissue. Fig. 6.10 shows this problem schematically. Fig. 6.10-A shows an initial state of the coupled system, in which a level-set formulation is used to divide the domain into the scaffold and medium parts (the $\psi$ function on the top) while a phase-field function divides the domain into tissue and void space (the $\phi$ function on the bottom). The positive side of the level-set function ( $\psi>0$ ) is not defined in the domain of the tissue growth model ( $\phi=\text{NaN}$ ). After several time steps, the system enters a state similar to the one depicted in Fig. 6.10-B, where the scaffold has shrunk to the left due to biodegradation, and the tissue has grown to the right. As can be seen in the figure, these movements cause a new area to appear between the scaffold and neotissue, in which the level-set function is negative ( $\psi<0$ ) and the phase-field variable is one ( $\phi=+1$ ). Since this area is part of the void space for the phase-field model, the tissue should grow in this direction as well, but this behavior does not have a clear definition from an implementation point of view. From a common-sense perspective, one might expect part of the tissue to still be attached to the surface of the scaffold and move with it to the left, a behavior that is quite challenging to consider in the formulation of the two coupled interface moving problems. Moreover, it is not fully clear how such behavior happens in experiments, so dedicated experiments are required to observe the actual reaction of the growing neotissue to the shrinkage of the scaffold. This behavior may be sensitive to experimental conditions such as static or perfusion setup of the medium, making the problem even more complex from the mathematical perspective.

Figure 6.10: Schematic presentation of coupled neotissue growth and biodegradation models (implemented using phase-field and level-set methods respectively): A) initial state of the system with level-set ( $\psi$ ) and phase-field ( $\phi$ ) variables dividing the domain for the biodegradation and neotissue growth models, respectively, B) the state of the system after some time steps, in which the biodegradation shrinks the scaffold to the left, the neotissue grows to the right, and a new area emerges between the scaffold and the neotissue.

A potential solution for the mentioned coupling challenge would be to take advantage of the multi-phase field method for describing both the biodegradation and tissue growth processes using a single model. Putting this solution into practice requires a tremendous amount of work to re-implement the biodegradation model using the phase-field method, which needs to change part of the fundamental equations in the mathematical model elaborated in Chapter 3, but it removes the necessity of coupling two models and dealing with the emerging difficulties in the implementation. Such an idea has already been implemented in studies like Moure et al. work [199, 200], in which they used a multi-phase field model to investigate individual and collective cell migration and crawling.

The next challenge would be the validation of the coupled model. In addition to the necessity mentioned above for performing dedicated experiments to observe the actual behavior of tissue growth on biodegradable scaffolds, such experiments are crucial to validate the model from a quantitative point of view. Doing these experimental studies seems to be challenging and resource-demanding because both the qualitative behavior and quantitative measurements should be recorded. The qualitative behavior is crucial to observe how tissue growth reacts to the shrinkage of the material underneath, while the quantitative output of the experiments can be directly used to validate the rate of degradation and neotissue formation predicted by the coupled model. Such resource-demanding experiments were never planned as part of this PhD, and as a result, the coupled models could not be easily validated.

A suitable workaround for the aforementioned challenges can be coupling the biodegradation model with a simpler tissue growth model in the first place. A simpler model here implies that it does not include an interface tracking sub-model, for which the above challenges may appear in the implementation. An example of such simplified coupling is the work of Byrne et al. [43], in which random walk algorithms were used to model neotissue differentiation, scaffold degradation, and their coupled effect phenomenologically. Another example of a more straightforward model is Carlier et al. work [44, 45], in which they solved a set of taxis-diffusion-reaction PDEs describing the evolution of biochemical factors, cells, and extracellular matrices to model the healing process, assuming the entire domain is already occupied by some form of tissue (no void space). Due to the lack of interface capturing, such a model is easier to integrate with the biodegradation model. This approach can be more difficult to validate, but it comes with the advantage of a more straightforward implementation. Looking at the time scales of the both processes however, several weeks for neotissue filling of a 3D porous scaffold of 6 mm cube and several months to years for the degradation of the same scaffold (depending on composition and circumstances), a practical work-around might be to decouple them. This would mean first simulating the neotissue growth until complete filling and only then considering scaffold degradation which focuses all the attention on the zone between the tissue and the scaffold.

6.9 Conclusion

In this chapter, two curvature-driven models of the neotissue growth process were developed using the phase-field and level-set interface tracking methods. The models were used in various 2D and 3D growth simulations with cells seeded on the circumference of basic geometrical primitives, showing the similarities and differences between the phase-field and level-set methods. Additionally, the level-set model was used to simulate tissue growth on open porous scaffolds generated with TPMS geometries, demonstrating the applicability of the developed models for tissue engineering applications. The complications to couple the tissue growth models with the developed biodegradation model in Chapter 3 were also elaborated, in which the potential solutions to address the emerging coupling challenges were discussed.

Chapter 6 Computational modeling of the neotissue growth process