## Abstract

Tissue engineering is becoming consolidated in the biomedical field as one of the most promising strategies in tissue repair and regenerative medicine. Within this discipline, bone tissue engineering involves the use of cell-loaded porous biomaterials, i.e. bioscaffolds, to promote bone tissue regeneration in bone defects or diseases such as osteoporosis, although it has not yet been incorporated into daily clinical practice. The overall success of a particular bone tissue engineering application depends strongly on scaffold design parameters, which do away with long and expensive clinical protocols. Computer simulation is a useful tool that may reduce animal experiments and help to identify optimal patient-specific designs after concise model validation. In this paper, we present a novel mathematical approach to bone regeneration within scaffolds, based on a multiscale framework. Results are presented over an actual scaffold microstructure, showing the potential of computer simulation, and how it can aid in the task of making bone tissue engineering a reality in clinical practice.

## 1. Introduction

The aim of bone tissue engineering, as a subdiscipline of tissue engineering, is to regenerate bone mass loss due to tumour or any other physiological disorder, especially in large bone defects, to accelerate bone fracture healing or to promote bone growth under certain circumstances or diseases such as osteoporosis. This subdiscipline has received an increasing interest since the early works presented in the 1970s (Levin *et al*. 1975; Rejda *et al*. 1977). In the last decade, a wide number of experimental works using scaffolds for bone regeneration have been presented. The strategy commonly involves the use of cells (usually bone marrow stromal cells (BMSC), mesenchymal stem cells (MSCs) or preosteoblasts), combined with bone morphogenetic proteins or growth factors such as transforming growth factor-β (TGF-β), seeded within porous scaffolds (usually made of polymeric, bioceramic or bioglass materials) to promote bone tissue after implantation. The scaffold's behaviour *in vivo* has been evaluated in many papers, although numerical simulation may help in the understanding of such behaviour with the use of the computer.

Very few works have been conducted within the framework of numerical simulation in tissue engineering. It should be stated that the main objective of mathematical modelling and numerical simulation in this field is to predict the behaviour of scaffolds in order to get an optimal design according to its functionality. In tissue engineering methodology (see Evans *et al*. 2007), scaffolds may be preseeded with cells *in vitro* prior to implantation in the human body. In this step, an appropriate pore size, pore interconnection, porosity and permeability must be assured in order to permit cells to move through the scaffold core and adhere onto its surface. In bone tissue engineering applications, it is important that the scaffold is strong enough to support the mechanical loads transmitted by muscles and joints, which will be brought about by its effective (macroscopic) mechanical properties. Both permeability and effective mechanical properties may be tailored individually by gaining control over some scaffold parameters microscopically, such as mean pore size, porosity and bulk biomaterial mechanical properties. In this context, Hollister and co-workers (Hollister *et al*. 2002; Taboas *et al*. 2003; Lin *et al*. 2004) focused on the design of bone scaffolds as an optimization problem to obtain a microstructure as similar as possible (mechanical properties, porosity, pore size, etc.) to that of the implanted region. For this purpose, homogenization theory was extensively applied in their studies. Additionally, mechanical properties and permeability were homogenized over a representative element of a bioceramic scaffold microstructure (Sanz-Herrera *et al*. 2008*c*). The obtained results were corroborated by an experimental set-up showing the potential of numerical tools for the characterization of scaffold properties.

Ideally, the scaffold should degrade, giving way over time to new bone regeneration and therefore complete scaffold disappearance. Scaffold degradation usually takes place by chemical pathways (hydrolysis) in the case of polymeric scaffolds. This issue has been the subject of numerical analysis by Adachi *et al*. (2006) and recently by Wang *et al*. (2008). A proper understanding of hydrolysis kinetics can modulate the rate of scaffold mass loss, allowing us to choose a specific polymer biomaterial for a specific application. Other issues related to scaffold microarchitecture, such as nutrient supply and waste removal, have been modelled computationally and may be reviewed in the work conducted by Sengers *et al*. (2007).

Cell–biomaterial interaction is considered to play a major role in tissue engineering. In recent years, it has been observed that adherent cells sense and respond to the mechanical environment they are anchored to (Discher *et al*. 2005; Engler *et al*. 2007). The mechanical environment is mainly related to substrate stiffness, although topological factors such as surface roughness or curvature are of extreme importance as well (Pelham & Wang 1997; Wong *et al*. 2003; Engler *et al*. 2004). Moreover, mechanical forces in the intra- and extracellular space can induce a wide range of cell processes such as growth, differentiation, apoptosis, migration, gene expression and signal transduction (Lauffenburger & Horwitz 1996; Kitano 2002; Bischofs & Schwarz 2003; Engler *et al*. 2007; Ghosh *et al*. 2007; Suresh 2007). In fact, cells transduce mechanical stimuli into biochemical signals and vice versa, i.e. mechanotransduction. Cell migration, proliferation and differentiation have a prominent significance in bone tissue engineering processes (Jerome & Maquet 1997; Ishaug-Riley *et al*. 1998) both *in vivo* and *in vitro* and are mainly affected by scaffold stiffness, microarchitecture and biomaterial composition. All these arguments have drawn the attention of mechanicians focusing their interests on another emerging research line, namely, *cell mechanics*. There is today a vast amount of literature specially relating to mathematical modelling and computer simulation (DiMilla *et al*. 1991; Kitano 2002; Bischofs & Schwarz 2003; Nicolas *et al*. 2004; Zaman *et al*. 2005; Maini *et al*. 2006; Schwarz *et al*. 2006; Suresh 2007; Moreo *et al*. 2008) showing significant results and direct application to tissue engineering.

From all the above, it becomes clear that bone tissue regeneration *in vivo* using scaffolds inherently attends to two well-differentiated spatial and temporal scales. One scale is at the tissue level or macroscopic scale, whereas the other is at the pore level or microscopic scale. Traditionally, mathematical modelling of tissue growth within scaffolds has been restricted to one of these scales. For instance, at the macroscopic scale, osteochondral defect repair using scaffolds has been simulated by Kelly & Prendergast (2006) using a previous mathematical model of cell/tissue differentiation (Prendergast *et al*. 1997). A macroscopic approach for the simulation of bone tissue regeneration within scaffolds considering the effect of its microstructure was introduced in our previous work by Sanz-Herrera *et al*. (2008*b*). At the microscopic scale, Byrne *et al*. (2007) presented a mechanobiological model applied to a periodic unit cell of the scaffold microstructure. In that model, the influence of several factors such as porosity, mechanical properties and others were analysed with regard to tissue regeneration. At the same microscale, bone tissue regeneration within a scaffold unit cell was numerically reproduced by Adachi *et al*. (2006) using concepts and hypotheses previously established for a bone remodelling theory (Adachi *et al*. 1997, 2001). Recently, a parametric analysis of several factors of major importance in tissue engineering has been presented by Sanz-Herrera *et al*. (2009). The model was based on a multiscale technique applied to a bone tissue engineering example. The underlying microstructure of the scaffold used in the model was simplified by considering a periodic face-centred cubic repetition of the pores in the unit cell (Brígido-Diego *et al*. 2007).

This paper focuses on the mathematical modelling and computer simulation of bone tissue engineering events in the presence of scaffolds (see table 1 for nomenclature). This work differs from the previous ones in that it employs a mathematical methodology that considers two scales of analysis. This is critical in bone tissue engineering since scaffold implantation occurs macroscopically, although the most important biophysical phenomena and control design parameters stand microscopically. In our model, the macroscopic level includes the bone organ and scaffold region, whereas the microscopic scale is associated to the pore scaffold level. Even when the presented model considers multiphysics phenomena, the main underlying hypothesis is that bone regeneration is regulated by mechanical factors.

The proposed multiscale model is described mathematically in §2. An example of application is analysed in §3, showing results at both the macroscopic and microscopic levels. The example includes a realistic microstructure of a polymeric porous biomaterial. Finally, the paper closes with some concluding remarks regarding the perspectives and future directions of the bone tissue engineering field from a clinical perspective, and how this work and computer simulation can provide qualitative information in parallel to animal experimentation.

## 2. Mathematical methodology

The mathematical model based on a multiscale methodology is presented in this section. We first introduce the governing equations and their weak form, distinguishing between bone organ and scaffold domain. The scaffold macro-domain is coupled with its underlying microstructure. Moreover, the internal microstructure of the scaffold evolves as cells form new bone tissue. The modelling of this phenomenon as well as the coupling between the micro- and macro-scales of the problem at the scaffold region is explained in §2*b*. Finally, the macroscopic model of the bone organ is introduced in §2*c*.

### (a) Multiscale analysis

Let us consider the macroscopic region to be composed of two subdomains, namely, the scaffold region and the non-scaffold region , i.e. the bone organ. The domain labelling may be seen in figure 1, which in turn will serve as application (§3) to the developed model. We assume that bone remodelling takes place at the bone organ region, while microscopic bone apposition occurs at the microsurface of the scaffold. As hypothesis, we consider that these phenomena are governed mainly by mechanics through a mechanical stimulus, similar to other bone remodelling theories (Pauwels 1965; Cowin & Hegedus 1976; Carter *et al*. 1987; Huiskes *et al*. 1987; Beaupré *et al*. 1990; Jacobs 1994). The macroscopic elasticity equations in the absence of body forces are posed both at the scaffold and bone organ as(2.1)with ** ϵ**(

**) and**

*x***(**

*σ***) being second-order tensor fields denoting the (linearized) strain and stress, respectively, and**

*x*^{0}the macroscopic fourth-order elasticity tensor defined as(2.2)where

**denotes the macroscopic scale and**

*x***the microscopic one. Note that the elasticity tensor is a function of each macroscopic point, being therefore a heterogeneous material. In fact, it will be a function of a field variable (apparent density) defined in the bone organ as explained below. However, the effective mechanical properties of the scaffold are directly derived from the evolution of the underlying microstructure. This evolution is determined by the new bone apposition, scaffold degradation and invasion of cells to the scaffold core. The model of microscopic bone growth within the scaffold and biomaterial degradation is explained below.**

*y*As pointed out, we macroscopically model the process of cell migration and to some extent vasculogenesis in the scaffold region. This is based on the assumption that bone regeneration follows the same sequence of events as conventional bone regeneration under normal biocompatibility conditions. First, a haematoma is formed and inflammation occurs (Einhorn *et al*. 1995). At this level, the walls of the scaffold are filled with MSCs and preosteoblastic cells. This migration process from the surrounding bone to the scaffold core is here modelled by a standard Fick's law,(2.3)where the diffusion matrix denoted by ** D**(

**,**

*x***) is considered fully anisotropic and ,**

*y**c*(

**) being the number of cells per unit volume in the scaffold domain and**

*x**c*

_{s}the number of cells per unit volume surrounding the scaffold boundary, which is assumed to be the maximum density of MSCs available. Note that, as with mechanical properties, the diffusion coefficient is obtained from an analysis of the scaffold microstructure, which couples together the macro- and micro-scales.

With (2.1) and (2.3) in hand, we can apply the variational form of these equations as(2.4)(2.5)the space of admissible variations being defined as(2.6)with *H*^{1}(*Ω*^{0}) being the first-order Sobolev space. The derivation of ** D** is explained below.

For simplification, we have considered that the bone organ and the scaffold are fully connected. Moreover, the free boundary of the scaffold (which is not connected with the rest of the bone organ) does not support boundary loads. Using a homogenization procedure in (2.5), we obtain the two-scale variational form of the mechanical problem, namely(2.7)with ^{0} obtained through the homogenization process as explained below.

Note that (2.7) is a direct consequence of the linearity of the problem. In nonlinear situations, the separation between the macro- and micro-scales cannot be accomplished unless an appropriate linearization is conducted (Terada & Kikuchi 2001). This statement is valid in classical multiscale linear mechanics, where no evolution of the microstructure occurs. Nevertheless, since the considered microstructure is expected to grow and to degrade, an additional hypothesis is introduced. We consider that both bone growth and scaffold degradation are averaged during each time iteration. This implicitly considers that the macroscopic apparent properties are unchanged during the time step. A consequence of this simplification is that the time scales associated with the application of loads and diffusive analyses are uncoupled from those of bone ingrowth and degradation.

Equations (2.4) and (2.7) are numerically implemented in a finite-element framework using the Abaqus software (HKS 2001). Parallel performance is used to accelerate the micro–macro approach. The details of the numerical implementation are described by Sanz-Herrera *et al*. (2008*a*).

### (b) Model for the scaffold domain

The scaffold domain corresponds to the underlying microstructure of a macroscopic material point . It is treated in the form of a representative volume element (RVE) that statistically represents the microstructure of such a macroscopic point. We distinguish between fluid and solid domains in the RVE, such that , with and being the microscopic solid and fluid domains, respectively (see figure 1). The microscopic scale is denoted by ** y**.

Bone ingrowth occurs at the scaffold surface, i.e. , after MSC attachment and subsequent proliferation, differentiation and bone matrix production. We assume that bone formation is driven by a mechanical stimulus evaluated at the microstructural level and given through macroscopic quantities. Furthermore, we consider that the activity of bone cells and subsequent bone matrix production and/or resorption happen at the scaffold pore surface. Then, we establish that the rate of bone mass formation at the local scaffold microsurface is(2.8)where the bone mass rate is a function of a microscopic mechanical stimulus *Ψ*^{φ} and the current macroscopic population of MSCs, . The derivation of the microscopic stimulus from macroscopic information is addressed below. Additionally, we consider that the activity of the cells attached to the scaffold surface starts when a mechanical threshold is achieved. Above this threshold, MSCs start to differentiate into osteoblasts and deposit bone matrix on the surface, while under that value they preferentially stay inactive (Turner *et al*. 1994; Tang *et al*. 2006). Therefore, we can express bone growth on the scaffold surface within the scaffold domain in the following way:(2.9)where is an empirical rate function (*k* being a model parameter) proportional to the macroscopic concentration of MSCs, , currently attached to the surface of the scaffold. The assumed microscopic mechanical stimulus *Ψ*^{φ} is analogous to that in the previously established theories of bone growth and remodelling (Beaupré *et al*. 1990), but conceptually different since we consider it here to be a microscopic stimulus evaluated at the scaffold pore surface. This is stated as(2.10)where *N* is the number of different load cases; *n*_{i} is the average number of cycles per time unit for the load case *i*; *m* is a model parameter (Beaupré *et al*. 1990); and is an effective scalar stress expressed as(2.11)with local Young's modulus *E* distinguishing among scaffold (s), non-mature bone (n) or mature bone (m). Finally, *U*^{φ}=1/2*σ*^{φ} : *ϵ*^{φ} is the strain energy density at the scaffold microsurface. This is obtained through the localization problem and subsequent distribution of strains along the scaffold microstructure described below.

As a first simplified approach, the time evolution of the mechanical properties of non-mature bone is considered to be linear (Adachi *et al*. 2006),(2.12)where , are the fourth-order elasticity tensors of non-mature (n) and mature bone (m), respectively. Note that, at the microscopic level, the mechanical properties depend only on the mineralization ratio since the distinction between cancellous and cortical bone does not make sense at this level. The newly formed bone matrix is considered at this scale as isotropic, as a first approach.

The degradation mechanism is modelled as a hydrolysis process within the solid microstructure of the scaffold. This model was presented by Adachi *et al*. (2006). Briefly, the water content in the polymer chemically reacts and produces an internal change of the polymer, causing it to detach, consequently leading to erosion of the biomaterial bulk (Gopferich 1997). We consider that the dimensionless spatial rate of water content *d*, defined as the ratio between the concentration of water (mol) in the bulk of the polymer and the concentration of water (mol) at the boundary of the scaffold, is governed by the diffusion equation (Adachi *et al*. 2006)(2.13)with *α*>0 being the diffusion coefficient. On the other hand, the rate of change of the molecular weight of the polymer due to hydrolysis is assumed to depend on the local water content (Adachi *et al*. 2006) (0≤*d*≤1) as(2.14)with *β*>0 being a material constant. The mechanical properties of the polymer are assumed to be linearly related to its molecular weight (Adachi *et al*. 2006) as(2.15)with *W*_{0} being the reference molecular weight of the polymer and being the associated elasticity tensor.

The bone growth process is numerically modelled by using the voxel finite-element method (FEM; Adachi *et al*. 2006). Thus, the amount of bone to be formed within the scaffold is simulated by adding voxels, representing the newly created bone (Sanz-Herrera *et al*. 2008*a*). As a simplification, we considered that bone formation occurs only on the surface of the scaffold. The mechanical stimulus in equation (2.10) is computed at the centroid of each voxel to be the candidate for bone formation, i.e. the boundary voxels. Then, equation (2.9) is integrated by the forward Euler method and the number of voxels corresponding to the new bone volume production *M*_{c}/*ρ*_{b} are added, where *ρ*_{b} is the mineral bone density and *M*_{c} a model parameter, during each period or time step. For the implementation, the voxel scaffold microstructure is treated as a binary third-order matrix with indices *i*, *j* and *k* (with *i*, *j*, *k* from 1 to the considered length of voxels). The matrix contains ‘1’ at the spatial position *ijk* when a bone or scaffold voxel is identified at that position, and ‘0’ when it corresponds to a void voxel. Then, we identify new voxels to be formed when a shift between 0 and 1 is computed along the scaffold microstructure. Then, the position that contains a 0 automatically changes to 1, meaning that a new bone voxel has been created. Analogously, scaffold degradation is simulated within a voxel FEM framework. Here, we remove the voxels where *W*<*W*_{c}, *W*_{c} being a model parameter, during each time step. A scheme of this procedure can be seen in figure 2.

The remainder of this subsection is devoted to the analysis of the RVE in order to obtain the homogenized properties ^{0} and ** D** available for the implementation of equations (2.4) and (2.7). Moreover, this analysis allows us to obtain the microscopic distributions of the field variables involved. The linear, elastic microscopic problem in the absence of body loads in the solid domain in reads as(2.16)with 〈.〉 being the volume average over the RVE. The absence of boundary conditions in (2.16) is noticeable. These boundary conditions must reproduce, as closely as possible, the

*in situ*state of the RVE inside the material. Therefore, they strongly depend on the choice of the RVE itself, and especially on its size. Two classical boundary conditions have been established traditionally,(2.17)with

**being the outward normal to the boundary. Equations (2.17) refer to uniform stresses and strains, respectively, on the boundaries of the solid domain of the RVE. It is immediately seen that a divergence-free field**

*n*

*σ*^{φ}that satisfies (2.17)

_{1}also satisfies 〈

*σ*^{φ}〉=

*σ*^{0}, and a displacement field that satisfies (2.17)

_{2}also satisfies 〈

*ϵ*^{φ}〉=

*ϵ*^{0}(see Suquet (1987) for further reading).

If periodic media are under consideration, the boundary conditions (2.17) must be rejected (Suquet 1987). Therefore, the periodicity conditions imply the following statements:

the stress vectors

*σ*^{φ}.are opposite on opposite sides of the contour , and*n*the local strain

*ϵ*^{φ}(*u*^{φ}) is split into its average and fluctuating terms such that(2.18)

Thus, the periodicity boundary conditions yield(2.19)

Let us substitute equations (2.18) into the equilibrium (divergence-free) equation of the stress tensor (see equation (2.16)_{1}), namely,(2.20)where the boundary conditions are related to (2.17) or (2.19). By virtue of the linearity of the problem, the solution of in equations (2.16) for a general macro-strain *ϵ*^{0} may be expressed as the superposition of elementary unit strain solutions *ϵ*^{φ}(*Χ*_{kh}) (Suquet 1987), such as(2.21)where *Χ*_{kh} are the displacements associated with those elementary strain states denoted by indices *kh* resulting from the solution of equations (2.16).

By substitution of equation (2.21) into equation (2.18)_{1}, the microstrains are expressed as(2.22)where *I*_{kh} is the fourth-order identity tensor with components . Equation (2.22) allows us to obtain the distribution of the strain field along the microstructural domain from the macroscopic strain field. For that purpose, proper boundary conditions have to be previously established according to (2.19). Then, we can obtain the strain energy density from the microscopic stress and strain fields to evaluate the micromechanical stimulus needed in equation (2.11) (Sanz-Herrera *et al*. 2008*a*).

Once the elementary solutions *ϵ*^{φ}(*Χ*_{kh}) through (2.22) have been obtained, the macroscopic stress–strain relationship is obtained straightforwardly as follows:(2.23)

Consequently, at the macroscopic scale, the elasticity tensor is identified in equation (2.23) as(2.24)

Using the standard formulation, the variational form of equation (2.20) yields (Suquet 1987)(2.25)where the space *V*_{Y} is defined as(2.26)with *H*^{1}(*Ω*^{φ}) being the first-order Sobolev space. Using (2.21), equation (2.25) can be further developed, namely(2.27)where represents the characteristic microstructure displacement in the *p*-direction due to an applied *kl* unit strain, *k*=*l* being the normal unit strain states and *k*≠*l* being the shear unit strain states; *ω*_{i} represents a virtual displacement. There are six strain states in total (three normal and three shear) corresponding to the six linear equations above (equation (2.27)). Once these functions have been obtained, the macroscopic stiffness tensor can be computed through equation (2.24).

Homogenization theory is used in fluids as well to obtain the overall response of a microscopic flow in a macroscopically porous medium, i.e. Darcy's law, and the microscopic distribution of the velocity field. Therefore, here we establish an analysis similar to that of solids to analyse domain . The equilibrium and continuity microscopic equations for an incompressible viscous fluid within the steady-state Stokes' flow problem read as(2.28)with *p*^{φ} being the fluid pressure; *v*^{φ} the fluid velocity; *ρ*^{φ} and *μ*^{φ} the fluid density and dynamic fluid viscosity, respectively. As in the solid domain, we consider the overall velocity to be the result of the macroscopic velocity field plus the perturbation induced by the microstructure, i.e. , with *φ*=*l*/*L* being the ratio between the micro and macro levels. The body forces vector per unit volume is denoted by ** f**. For brevity, we define , with

**1**being the second-order identity tensor. Since we are dealing with periodic boundary conditions, we establish an analogous procedure to (2.19) for the fluid problem,(2.29)

As in the solid problem, the microvelocities in (2.28) admit a solution of the form(2.30)with being the characteristic fluid velocities associated with unit pressure gradients (Terada *et al*. 1998). Equation (2.30) permits us to obtain the microvelocities along the microscopic domain through the macroscopic pressure.

After some algebraic manipulations, and using (2.28) and (2.30), functions (Terada *et al*. 1998) must obey the following expression in the microscopic fluid domain of the RVE:(2.31)in combination with the boundary conditions (2.29).

To obtain the homogenized macroscopic Darcy's law, the microvelocities of (2.30) are averaged over the microstructural domain (RVE) as follows:(2.32)with ** K** being the macroscopic permeability matrix. To derive (2.32), the null fluctuation of the velocity field is taken into consideration. The Fick's diffusion macroscopic matrix is related to the permeability as

**=**

*D***.**

*K**D*

_{0},

*D*

_{0}being a model parameter.

The variational form of the fluid phase is developed analogously to that of solids, as shown before. Therefore, by using periodicity assumptions in (2.29), the variational form of equations (2.28) yields(2.33)with the space *V*_{Y} being defined as earlier in this section. Using the superposition of (2.30) and the result (2.31), we finally get the following integral expression:(2.34)where represents the characteristic microstructure fluid velocity due to an applied generalized pressure gradient in the direction *i* and where *γ*_{j} represents a virtual velocity. Once these functions have been obtained, the macroscopic permeability *K*_{ij} may be obtained from equation (2.32).

### (c) Model for the bone organ domain

The bone organ is only treated in macroscopic terms by simulating the phenomenon of bone remodelling. Carter's group (Carter *et al*. 1989; Beaupré *et al*. 1990; Jacobs 1994) proposed an isotropic remodelling model with the mechanical stimulus identified with the so-called daily tissue stress level, *Ψ*^{0}, a scalar quantity similar to the one defined by Carter (Carter *et al*. 1987) but generalized to include several load cases. It is defined as(2.35)where *N* is the number of different load cases; *n*_{i} the average number of load cycles per time step for each load case *i*; *m* an experimental parameter varying from 3 to 8 (Whalen & Carter 1988); and the so-called effective stress at the tissue level.

The first assumption of this model is that the order of load application has no influence on the adaptive response of the bone. This seems reasonable owing to the difference between the time scales of causes (load frequencies) and effects (remodelling process) as computationally proved by Jacobs (1994). Beaupré *et al*. (1990) also established a relation between the stress at the continuum level and the one at the tissue level , following a standard homogenization procedure and supported by experimental data. With these two assumptions, *Ψ*^{0} may be rewritten as a function of stress at the continuum level as(2.36)with being defined by(2.37)where *E* is the elasticity modulus and *U*^{0} the internal strain energy at the macroscopic level; *ρ* is the apparent density; and is the density of the ideal bone tissue with null porosity.

The next step corresponds to the establishment of the evolution law for density as a function of the external stimulus. Bone resorption and apposition occur at internal surfaces, so the mechanical stimulus is directly related to the surface remodelling rate, , which quantifies the bone volume formed or eliminated over the available surface per time unit. Following this idea, and as an approximation to the experimental curves, Jacobs (1994) proposed the following expression:(2.38)*k*_{i} being the slopes of the two straight lines of apposition and resorption, respectively; *Ψ*^{0*} a reference value of the tissue stress level; and 2*w* the width of the so-called ‘dead’ or ‘lazy’ zone. Consequently, the evolution law for the apparent density is expressed as(2.39)where the added or removed bone is assumed to be completely mineralized, that is, with a maximum density ; *S*_{v} is the specific surface; and is the ratio between the available surface for remodelling and the total internal surface (Jacobs' model considers that the whole existing surface is active, that is, , which leads to an excessively fast remodelling process). The elasticity modulus and Poisson's ratio are finally obtained by the following experimental correlation (Jacobs 1994):(2.40)The parameters we used regarding this model were obtained from Jacobs (1994).

## 3. Results

The mathematical model described above was mathematically implemented for an application example found in the literature (Savarino *et al*. 2007). A healthy left femur of a New Zealand rabbit was modelled by computer tomography (CT) using the Mimics software (Materialise 2006; figure 1*a*). On this solid model, a virtual segmental 3×7 mm defect in the proximal third of the left femur, laterally to the greater trochanter (Savarino *et al*. 2007), was performed in order to simulate the scaffold implantation. The FE mesh of that model was obtained using Harpoon software (Sharc 2006) with four-node bilinear tetrahedra in both the bone and scaffold macroscopic domains. Femoral loads were applied at the proximal part considering the rabbit hopping according to Gushue *et al*. (2005). Moreover, the initial femur density was obtained by applying the remodelling model presented above (Beaupré *et al*. 1990) in the healthy femur until convergence.

Since the mathematical model accounts for two different scales of analysis, the scaffold microstructure needed to be simulated as well. Hence, the actual geometry of a poly-*ϵ*-caprolactone scaffold was obtained through micro-CT. An RVE of such microstructure was selected, distinguishing between the solid (figure 1*b*) and fluid (figure 1*c*) parts. For simplification, we considered that this RVE represents the microstructure of the entire biomaterial, such that its microstructure may be reproduced by a periodic repetition of such RVE. The scaffold microstructures were meshed using voxel FEs for the purpose of simplicity (Adachi *et al*. 2006; Sanz-Herrera *et al*. 2008*a*). The model parameters used for this simulation are shown in table 2, whereas the scaffold properties are highlighted in table 3.

The results of the simulation of the example described above, obtained by implementing the multi-scale model, are represented in figure 3. We can observe the apparent bone density distribution, the macroscopic results in the bone region in terms of apparent density and the microscopic new bone deposition distribution along the scaffold architecture. The microscopic results corresponding to any point could be illustrated. For example, the results at the macroscopic midpoint of the scaffold biomaterial are shown in figure 3. This figure also shows bone mass evolution within the scaffold microstructure over time and the resorption of the scaffold in the last steps of the analysis. Note that at the microscopic level new bone formation and apposition occur at surfaces where a higher mechanical stimulus is found, according to our model. Since this variable depends on scaffold microstructure, it can be concluded that the scaffold microarchitecture is key to new bone tissue formation. Within the context of the influence of scaffold morphology on bone regeneration, further information can be found in the experimental studies carried out by Tsuruga *et al*. (1997), Kuboki *et al*. (2001) and Karageorgiou & Kaplan (2005), or as recently presented by Andreykiv *et al*. (2008) within the framework of the effect of peri-implant surface geometry, from a computational perspective. One of the conclusions derived from the analysis is that the effect of the scaffold implantation on the bone organ (in terms of bone remodelling) is negligible, except for a very localized area surrounding the scaffold proximity region.

The overall bone regeneration and scaffold resorption are depicted in figure 4. The former is computed as the percentage of new bone regeneration in the microstructure (averaged per volume unit) at each macroscopic point. The scaffold resorption is computed analogously but considering the biomaterial scaffold mass loss in this case. Figure 4 shows that bone regeneration starts around day 22, which means that the critical mass *M*_{c} has been reached after that implantation time, after which an increasing trend over time is obtained. This is connected with many experimental works where bone regeneration increases with implantation time. Owing to the dispersion and variety of applications, we should only aim at validating this model qualitatively under similar conditions. For example, Pothuaud *et al*. (2005) implanted 75 per cent of porous hydroxyapatite–tricalcium phosphate scaffolds in the proportion of 4 to 6 loaded with BMSC, in the femoral condyle of New Zealand rabbits. One group of the scaffolds used in this experiment were functionalized with a cyclo-RGD (arginine–glycine–aspartic acid) peptide to enhance cell adhesion onto the scaffold surface. Results showed 12.3 and 21.3 per cent overall bone regeneration for non-functionalized and functionalized scaffolds at two weeks, respectively; and 21.1 and 18.1 per cent overall bone regeneration for non-functionalized and functionalized scaffolds at four weeks. On the other hand, injectable calcium phosphate into the knee joint of New Zealand rabbits was presented by Gauthier *et al*. (2005). New bone apposition was evaluated through three-dimensional micro-tomographic imaging, resulting in approximately 30 per cent of the overall bone volume formation at six weeks after implantation. The results yielded by these experimental works are in the range of those reported in this work, as shown in figure 4. However, the microarchitecture, biomaterial composition and exact location of the implantation of the scaffold simulated in this work are different from those analysed in the previously mentioned experimental reports, which may explain the differences between the results. The scaffold resorption shown in figure 4 is quite abrupt, starting approximately at day 54, once the critical minimum molecular weight was encountered at the scaffold microstructure. This trend is typical in the degradation of polymeric materials and lies within the range of the experimental results presented by Sung *et al*. (2004). Nevertheless, the degradation mechanism of polymeric biomaterials, i.e. hydrolysis, is highly dependent on parameters regarding the internal scaffold microarchitecture such as pore size or porosity. Based on clinical evidence, the degradation rate of polymeric materials should be as slow as possible (Hutmacher 2000).

## 4. Concluding remarks

The use of numerical techniques in tissue engineering is a promising tool to reach clinical application and to help in the task of making tissue engineering a clinically viable reality. It allows useful information to be obtained in a large number of clinical cases with the use of the computer. This is key to tissue engineering, where the role of many factors and parameters are misunderstood or need to be elucidated. A deep understanding of each biological phenomenon has to be gained by modellers, and a strong collaboration of the fields involved is required.

In this scenario, we presented in this paper a novel mathematical approach available for modelling bone tissue regeneration using scaffolds. We also included the simulation of some realistic experiments of scaffolds implanted in rabbits. From these examples, we see that bone formation occurs gradually over time, whereas scaffold resorption starts abruptly, which is consistent with some previous experimental works (Sung *et al*. 2004; Gauthier *et al*. 2005; Pothuaud *et al*. 2005). The major implicit hypothesis of this model is that bone regeneration within the biomaterial is mainly driven by mechanical factors. This hypothesis is in agreement with the experimental evidence of bone regeneration in the presence of biomaterials found elsewhere (Bostrom *et al*. 1998; Simmons *et al*. 2001; van der Donk *et al*. 2002; Gardner *et al*. 2006). Furthermore, mechanical factors have also been considered to play an important role in bone regeneration within scaffolds, as shown in other theoretical models (Kelly & Prendergast 2006; Byrne *et al*. 2007). However, non-mechanical factors such as nutrients/oxygen supply or the scaffold's chemical surface, among others, may also be of extreme importance in these applications. In order to account for some of the most important phenomena, the model includes multiphysics as well as multiscale arguments for the proper simulation of the main biological processes. For example, it is known that, when a scaffold is implanted into a bone defect, it produces a haematoma in the area surrounding the scaffold region that contains osteoprogenitor and mesenchymal cells. Here, effects such as cell migration and vasculogenesis take place, which were modelled as a Fickean diffusion, as a first approach. A more concise analysis could be performed to better simulate this kind of biological event.

The need for including two scales of analysis is due to the fact that tissue regeneration and biomaterial degradation take place at the pore level, whereas the overall behaviour of the scaffold implant is better observed at the organ level, where bone–scaffold interaction occurs, and each scale determines the response of the other. Moreover, this approach allows us to properly model both the actual bone organ geometry and the scaffold microstructure, as demonstrated in the present paper in contrast to other works. Indeed, we only define two model parameters to be determined from experimental fitting. However, when actual scaffold patterns are considered in the simulation, as presented in this work, it is difficult to assume periodicity of the RVE, and periodic boundary conditions may be debatable. Nevertheless, even when the RVE is non-symmetric, periodic boundary conditions provide an overall good estimation for a proper RVE size (Suquet 1987). Moreover, the mathematical multiscale implementation results showed an unusual class of nonlinear macroscopic response, in the sense that variable fields evolve over time for constant boundary conditions. This fact comes from the evolution of the microstructure because of regeneration and degradation phenomena, this nonlinearity being unusual in classical multiscale methods applied to mechanics. Nevertheless, we only modelled the net amount of regenerated bone, which is a rough approximation of reality. Experimentally, cells form a biolayer onto the scaffold microsurface as they attach. After this, they differentiate into osteoblasts that fabricate primary bone, finally resulting in mineralized mature bone. A more complete modelling should include all these factors. Another point that may be improved is a finer geometrical approximation of bone growth and hydrolytic resorption as presented here by the voxel FEM.

The multiscale model requires further experimental validation, even though a great effort has been made in this work to find experimental applications in the literature. There are many input parameters of the model presented in this paper that are difficult to obtain from the literature, for instance the internal microstructure of the scaffold. This makes a simulation versus experiment comparison extremely difficult. An exhaustive validation of the proposed models should include multidisciplinary collaboration among physicians, veterinary surgeons, material and mechanical engineers as well as modellers. In this framework, we could get access to the main parameters of the experimental set-up, hence validation may be properly conducted. Once this task is performed, this work may be used for scaffold design in a specific patient. In this sense, the anthropometric and metabolic features of individual patients may be taken into account, as long as the results provided are only treated qualitatively.

Finally, it may be stated that tissue engineering is still a growing young field. It has not yet accomplished relevant or satisfactory results to be translated to clinical application. Many biomaterials need to be proven, protocols must be further tested and refined, and commercial interests have to prove this field to be economically viable. Nevertheless, given the scientific promise, potential social impact and the young age of the field, many believe that it should be only a matter of time until tissue engineering reaches the main stream of surgical practice (Fauza 2003).

## Acknowledgments

This work has been supported by the Ministerio de Educación y Ciencia of Spain (Project DPI2007-65601-C03-01). Their financial support is gratefully acknowledged. The reader may view our webpage at: http://i3a.unizar.es/gemm. There, one can find plenty of information about our group's activities and the models we are developing. Furthermore, the reader is kindly invited to contact J.M.G.-A. whenever necessary, to obtain additional information about the presented models.

## Footnotes

One contribution of 15 to a Theme Issue ‘The virtual physiological human: tools and applications I’.

- © 2009 The Royal Society