Abstract

This work introduces a computational method for designing bone scaffolds for maximum bone growth. A mechanobiological model of bone adaptation is used to compute the bone growth, taking into account the shape of the defect, the applied loading, and the existing density distribution of the bone in which the scaffold has been implanted. Numerical homogenization and a geometry projection technique are used to efficiently obtain surrogates of the effective elastic and diffusive properties of the scaffold as a function of the scaffold design and the bone density. These property surrogates are in turn used to perform bone adaptation simulations of the scaffold–bone system for a sampling of scaffold designs. Surrogates of the bone growth in the scaffold at the end of the simulated time and of the strain energy of the scaffold at implantation time are subsequently constructed from these simulations. Using these surrogates, we optimize the design of a scaffold implanted in a rabbit femur to maximize volume bone growth into the scaffold while ensuring a minimum stiffness at implantation. The results of the optimization demonstrate the effectiveness of the proposed method by showing that maximizing bone growth with a constraint on structural compliance renders scaffold designs with better bone growth than what would be obtained by only minimizing compliance.

1 Introduction

Bone scaffolds are porous materials designed to offer an alternative to natural bone for use in bone grafting [1]. Bone grafting is a surgical procedure in which donor bone is implanted to replace bone that is too severely compromised to heal on its own [2]. This damage may be caused by trauma or disease, or may be congenital. No matter the cause, the damage can be very debilitating [2,3]. The bone graft fills the damaged region or the bone defect. The graft then provides the necessary support for the growth of new bone in the defect [4]. Bone grafts are one of the most common transplantation procedures performed in the world, second only to blood transfusions globally [5]. More than four million procedures are performed annually worldwide, and demand is forecasted to increase as the populations of many industrialized nations age [6].

In the face of increasing demand, the current sources of graft material have serious limitations. The preferred source is the patient, in which case the procedure is known as an autograft [4]. However, the volume of bone available for autografts is inherently limited. Only so much material can be removed before the functionality of the donor bone is compromised, limiting their use to smaller defects [1,7]. Moreover, autografts require that the patient undergo a surgery to remove the bone in addition to the implantation surgery, adding the potential for complications at the donor site [1,7,8]. The primary alternative to autografts are allografts. In an allograft, the source of bone is a donor other than the person receiving the graft [7]. Allografts only require a single surgery on the patient and are not subject to the same defect size restrictions as autografts. However, allografts are plagued by a shortage of donor material, high costs, and elevated risks of disease transfer and rejection by the patient [2,7].

In light of the limitations of these two options, significant research efforts have been devoted to bone substitutes in recent years. Their clinical application remains limited, owing principally to difficulties in their design and manufacture [9,10]. This difficulty stems from the complex and frequently conflicting requirements any bone substitute must meet. Bone growth is a relatively demanding metabolic process. As such, the growth of new blood vessels into the defect is essential to facilitate the transport of nutrients, cells, and waste to and from the healthy bone surrounding the defect (the growth front) [11,12]. Thus, to ensure that the defect fully heals, a certain level of porosity must be maintained by the scaffold throughout the healing process. Simultaneously, a degree of mechanical compliance must be maintained, as bone growth typically requires stimulation in the form of mechanical loading to occur [11,13]. Scaffolds which are too stiff shield the bone from mechanical stress, inhibiting growth within the scaffold and potentially causing resorption of the surrounding bone [10,14]. Opposing these goals is the need for structural stability of the scaffold. Mechanical failure of the scaffold in vivo is not acceptable. Furthermore, excessive deformation impedes angiogenesis, slows the formation of bone, and results in a prolonged healing period [2,15]. Increasing the porosity of the scaffold for better angiogenesis and increased mechanical stimulation necessarily reduces its load bearing capacity and stiffness. Further complicating matters, every defect is different in terms of its shape and the loading to which it is subjected.

Previous works have used computational optimization techniques to design the scaffold and attempt to address these challenges. Efforts have proceeded along two main veins: the design of the unit cell of a periodic scaffold and the design of complete scaffolds. Works in the former vein typically assume the scaffold to be composed of a single unit cell that is repeated to fill up the defect region and attempt to find the optimal design for that unit cell (cf. Refs. [14,1621]). Works in the latter vein have sought instead to optimize the entire scaffold for a particular defect (cf. Refs. [2225]). This requires a more involved analysis but allows the structure to be optimized for the local stress distribution unique to a given anatomical location and defect shape.

Nearly all of these works have focused on the mechanical properties of the scaffold. The biological requirements of bone growth tend to be treated as structural constraints, if considered at all [9]. This approach overlooks some of the most important aspects of the scaffold’s design, since bone growth requires certain biological prerequisites to occur at all [11]. Moreover, even if the initial conditions of the scaffold are favorable for bone growth, that does not guarantee that they will remain so as the bone adapts [9]. Several authors have employed finite element-based mechanobiological simulations to investigate the effect of scaffold geometry on bone growth over time [2630]. These studies all show a strong link between the scaffold’s architecture and how well the bone heals. This link has been corroborated by animal studies [27,29].

In spite of the importance of considering the biological aspects of bone growth, the application of these models to scaffold design has been very limited to date. In Ref. [12], a scaffold unit cell is optimized so that the rate of bone growth is matched to the rate of scaffold resorption. The goal is to ensure an adequate measure of mechanical stimulus, and support is available throughout the bone growth process. Mechanobiological considerations are incorporated in the design of optimized scaffolds in Ref. [19], which focuses on the optimization of the rod diameter of a rhombicuboctahedral unit cell with the objective of maximizing bone growth before the combined bone–scaffold is so stiff that it inhibits further growth.

The focus of this work is to employ a full transient mechanobiological simulation of bone growth to optimize a ceramic scaffold composed of a periodic rectilinear lattice [31] (Fig. 1). The use of such a simulation allows the changes in mechanical and diffusivity properties of the scaffold resulting from bone growth over the course of the recovery period to be accounted for. The objective of the optimization is to maximize the bone growth volume fraction (defined as the fraction of the void space within the scaffold that is filled with bone) at the end of a specified growth period. This objective function is selected since it is directly related to the purpose of the scaffold, which is to support bone growth and consequently improve healing. A constraint on maximum structural compliance at the time of implantation is imposed to ensure that the scaffold has sufficient stiffness to carry the loads in the absence of bone growth. The design variables correspond to the geometric configuration of the scaffold, which can be controlled during manufacturing and implantation.

Fig. 1
An 8 mm diameter ceramic scaffold manufactured via DIW [32]
Fig. 1
An 8 mm diameter ceramic scaffold manufactured via DIW [32]
Close modal

To analyze the bone and scaffold, we obtain homogenized mechanical and diffusion properties of the lattice, for which we build a surrogate model as a function of the design variables. Using homogenized properties circumvents the need to generate a body-fitted mesh for each scaffold design. This is desirable because a body-fitted mesh with adequate element size would have a very large number of elements and would thus incur a large computational cost for the analysis. The construction of the properties surrogates is motivated by the computational cost of the homogenization procedure itself, since properties must be updated for every element in the scaffold mesh for every time-step of the transient simulation and for every iteration of the optimization. These surrogates can be inexpensively evaluated and allow for a highly efficient simulation. Moreover, they only need to be computed once for a given design configuration of the scaffold unit cell. Finally, a surrogate of the bone growth volume fraction as a function of the scaffold design parameters is created by performing the mechanobiological simulation for a sampling of scaffold designs. This surrogate can be exercised to understand the effects of the scaffold design parameters on the bone growth and the scaffold stiffness and to optimize the scaffold design. Using this surrogate model, the optimal scaffold parameters are then found and validated using the simulation.

The rest of this manuscript is organized as follows. Sections 2 and 3 present the details of the formulation and the computational implementation, respectively, including the construction of the surrogate of the effective properties of the scaffold, the construction of the surrogates for bone growth and scaffold stiffness, and the scaffold optimization. An example of the proposed methodology is presented in Sec. 4, and we draw conclusions of this work in Sec. 5.

2 Formulation and Implementation

Our approach to the periodic scaffold optimization is to use a two-level surrogate model as outlined in Sec. 1. The first surrogate model is of the effective properties of the scaffold, which we obtain via numerical homogenization (cf. Refs. [3335]). Instead of using a body-fitted mesh that conforms to the geometry of the scaffold rods, we use the geometry projection method [3638] to perform the analysis on a fixed grid. The geometry projection method allows us to move the rods and change their size freely without re-meshing [39]. This speeds the homogenization process and allows it to be fully automated. We detail this methodology in Sec. 2.3.

The second surrogate model is of the strain energy of the scaffold at the time of implantation and of the percent bone growth into the scaffold. The latter is obtained using a finite element-based mechanobiological simulation of bone growth in the scaffold. We adapt the bone growth model from the literature, which simulates bone growth in an implanted scaffold using a homogenized model of the elastic and diffusion properties of the scaffold. As mentioned in Sec. 1, methodology is computationally efficient because it avoids the need to use a conforming mesh of the scaffold in the simulation, which would lead to a large number of elements and consequently a substantial computational expense. The details of this second surrogate are provided in Sec. 2.4.

For both the homogenization and simulations, several common assumptions are made. We assume that both the bone and the scaffold material behave as linear isotropic materials in terms of diffusivity and elasticity; that no plastic deformation or brittle failure of the bone or scaffold occurs during the growth period; that the scaffold is bonded to the bone (i.e., no slippage or separation occurs at the scaffold–bone interface); that during the simulated healing period, both bone and scaffold resorption are negligible; that the diffusivity of bone is directly proportional to its normalized permeability [27]; and finally, that the diffusivity of the ceramic is negligible relative to that of the bone.

2.1 Scaffold Geometry.

We start our formulation by describing the geometric parameterization of the simple rectilinear scaffold considered in this work (Fig. 1). The rods in each layer are perpendicular to those of the adjacent layers. This gives the scaffolds a distinctive grid pattern if viewed perpendicular to the plane of the layers. This type of scaffold can be readily manufactured with biocompatible ceramic inks using robotic deposition techniques such as direct ink writing (DIW) [32]. In this work, we assume that the scaffold is constructed using such a robotic deposition method. Fundamentally, this method works by extruding (depositing) a colloidal ceramic suspension, or ink, out of a robotically positioned nozzle [40]. While extrudable, the suspension maintains its shape and can bridge small gaps without the aid of support material. The part structure is then built up layer by layer, similar to conventional fusion deposition modeling (FDM) printers for plastics. Once the structure is complete, the part is sintered. For our purposes, we assume that the distortion of the rod cross sections at rod intersections is negligible. Additionally, we restrict our consideration to scaffolds with uniform rod spacing and with rods that have a circular cross section. A unit cell of this type of scaffold is shown in Fig. 2.

Fig. 2
Unit cell of the rectilinear scaffold (enclosed by dashed lines) with rods shown in gray. z is the direction perpendicular to the deposited layers, and θ is the orientation about the z axis. The in-plane cross sections ((a) and (b)) are perpendicular to the xy plane and to each other. The out-of-plane cross section ((c) and (d)) is normal to the z axis. θ also corresponds to the angle between the axes of the rods in alternating layers and the x and y axes, respectively.
Fig. 2
Unit cell of the rectilinear scaffold (enclosed by dashed lines) with rods shown in gray. z is the direction perpendicular to the deposited layers, and θ is the orientation about the z axis. The in-plane cross sections ((a) and (b)) are perpendicular to the xy plane and to each other. The out-of-plane cross section ((c) and (d)) is normal to the z axis. θ also corresponds to the angle between the axes of the rods in alternating layers and the x and y axes, respectively.
Close modal
As in Ref. [31], two dimensionless values parameterize the geometry of the unit cell: the diameter-to-separation ratio
d/l:=dl
(1)
and the overlap fraction
α:=1ad
(2)
where d is the rod diameter, l is the in-plane separation between consecutive rods in the same layer, and a is the out-of-plane separation of rods in adjacent layers. Since the scaffold is sintered after deposition, we note that these overlapping regions bond the rods together. These bonds can be assumed fully effective in terms of mechanical properties. The third geometric parameter considered is the pore size ps. This parameter has units of length and gives the scaffold microstructure a scale relative to the overall size of the scaffold.

In addition to the three preceding geometric parameters, the orientation θ of the scaffold relative to the defect is also considered. For simplicity, only rotation about the out-of-plane axis is considered. Rectilinear scaffolds are stiffest along the axes of their rods, much in the same way that the lamina of a fiber reinforced composite is stiffest along the axis of the fibers. The stiffness rapidly decreases as the direction of loading rotates away from this axis, since rods are increasingly loaded in bending as opposed to tension or compression. The orientation parameter allows the effective stiffness of the scaffold to be tailored to the loading, even if the scaffold microstructure parameters are fixed, as shown in Ref. [22]. Examples of scaffolds with different values of the four geometric parameters are shown in Fig. 3.

Fig. 3
Rectilinear scaffold designs with different values of the geometric parameters. All the design samples in a row are of the same overall size and are identical except for the value of the variable explicitly noted. Units are not provided for ps as the purpose of this figure is only to provide a visual comparison of the effect of varying one design parameter at a time.
Fig. 3
Rectilinear scaffold designs with different values of the geometric parameters. All the design samples in a row are of the same overall size and are identical except for the value of the variable explicitly noted. Units are not provided for ps as the purpose of this figure is only to provide a visual comparison of the effect of varying one design parameter at a time.
Close modal

2.2 Bone Growth Model.

As previously mentioned, bone requires a certain amount of load or stimulus as well as an adequate supply of osteoblasts to grow. Osteoblasts are the cells which specialize in the repair and maintenance of bone, which are responsible for the actual growth process [11]. However, as new bone is deposited, the distribution of stresses in the bone and ease of fluid transport (i.e., diffusivity) change even under a constant stimulus. This means modeling the transient bone growth is required to characterize the performance of a scaffold. As mentioned in Sec. 1, the bone growth model adopted in this work is based on that of Refs. [11,27,41]. Some aspects of this model, as we detail in subsequent sections, are adapted to the rectilinear scaffolds considered here. This model uses strain energy density as a measure of mechanical stimulation, and diffusion as a model of cell transport. Based on the mechanical stimulus received and the relative concentration of cells at a point, the bone deposited there is estimated.

It is worth noting that there are other bone growth models which have been proposed, such as a curvature-based model (cf. Ref. [30]), a damage response-based model (cf. Ref. [42]), and a tissue-differentiation model (cf. Ref. [26]). These all take into account different factors in determining how bone grows. The curvature-based model determines the soft tissue growth rate on a surface based on its curvature and extrapolates bone growth from that [30]. The damage response model estimates the damage that occurs to bone under loading and uses that to estimate growth. This model allows bone to grow with anisotropic properties according to its stress state and allows for resorption of bone which is under-stressed [42]. The tissue differentiation model seeks to estimate how cells within a region will diffuse and differentiate, then calculates how much and what kind of tissue grows based on the cell type and concentration [26]. However, unlike all of these models, the one we use has been quantitatively verified in terms of bone growth [27].

The bone growth model consists of a finite element-based simulation that estimates the density change in bone on an element-wise basis. The bone is assumed to have an infinite supply of osteoblasts that is not depleted by the bone growth process. A linear static analysis on the bone–scaffold assembly is performed to determine the stress distribution, and a transient mass diffusion analysis is conducted to determine the concentration of osteoblasts in each element. At every time-step during the healing period, the density and material properties of each element are updated.

The analysis region Ω is separated into two domains: scaffold (Ωs) and bone (Ωb), with Ω=ΩsΩb and ΩsΩb=. A simple diagram of the domains is presented in Fig. 4. The bone domain contains all the pre-existing bones surrounding the scaffold. The scaffold domain consists of the scaffold itself and any bone that grows within it. The material properties of the two domains are calculated differently. Any elements in the bone domain are henceforth referred to as bone elements and any in the scaffold domain as scaffold elements. The calculations pertaining to these domains are covered in Sec. 2.2.1 for bone and Sec. 2.2.2 for the scaffold. However, the density update, which we now describe, is the same for both regions.

Fig. 4
Schematic of the bone and scaffold domains. The geometry of the bone and scaffold domains is arbitrary.
Fig. 4
Schematic of the bone and scaffold domains. The geometry of the bone and scaffold domains is arbitrary.
Close modal
As the density update is done on each element, we compute all relevant local quantities at the element centroid and assume that they are uniform within the element. Therefore, while the expressions shown in the following can be computed at any point x ∈ Ω, in the implementation they are computed at each element centroid. The density for each time-step is computed using the update
ρt+Δt(x)=min(ρ^b,ρt(x)+ρ˙(x)Δt)
(3)
where ρtt and ρt denote the bone density at times t + Δt and t, respectively, and ρ^b is the maximum possible density of bone, cf. Table 1.
Table 1

Bone growth model constants

VariableDescriptionValuesUnitsSource
mEmpirical weighting factor for mechanical stimulus4[11,46,47]
Ψ*Reference mechanical stimulus50MPa/6 h[27]
wDead zone half-width of mechanical stimulus12.5MPa/6 h[27]
csSensitivity to mechanical stimulus0.02μm/MPa[41]
ρ^bMaximum achievable bone density1.92g/cm3[41]
ρ~bMinimum allowed bone density0.05g/cm3[41]
SeffbEffective specific surface area for a bone element0.2[27]
SeffsEffective specific surface area of a scaffold element0.6[27]
D^bMaximum diffusivity of bone100mm2/6 h[27]
CKozeny constant of bone0.022[27]
D0Lower bound diffusivity0.001mm2/6 h
ξSmooth max parameter to avoid division by zero0.0001
VariableDescriptionValuesUnitsSource
mEmpirical weighting factor for mechanical stimulus4[11,46,47]
Ψ*Reference mechanical stimulus50MPa/6 h[27]
wDead zone half-width of mechanical stimulus12.5MPa/6 h[27]
csSensitivity to mechanical stimulus0.02μm/MPa[41]
ρ^bMaximum achievable bone density1.92g/cm3[41]
ρ~bMinimum allowed bone density0.05g/cm3[41]
SeffbEffective specific surface area for a bone element0.2[27]
SeffsEffective specific surface area of a scaffold element0.6[27]
D^bMaximum diffusivity of bone100mm2/6 h[27]
CKozeny constant of bone0.022[27]
D0Lower bound diffusivity0.001mm2/6 h
ξSmooth max parameter to avoid division by zero0.0001

It should be noted that the upper limit on bone density in (3) is an artifact of the numerical bone adaptation scheme. In vivo, the maximum bone density would be naturally limited by the depletion of the available surface for bone to grow on, and thus there would be no need to impose such a limit. In our numerical scheme, however, the approximation of the specific surface S available for bone to grow (discussed in Sec. 2.2.2) is not exactly zero when the bone density attains its maximum value ρ^b, and consequently the bone growth rate of (4) would not be exactly zero as it would be in vivo. Moreover, it is numerically possible that the magnitude of the bone density rate of (4) and the time-step Δt in (3) are such that the density update of (3) exceeds the maximum density value. To prevent these situations, in a similar manner as in Ref. [27], we impose the maximum density limit in (3).

A small minimum allowable bone density ρ~b is also imposed everywhere to ensure that the mechanical properties of bone are not zero, which leads to an ill-posed mechanical analysis. This situation can arise when determining the initial bone density from a micro-computed tomography (micro-CT) scan based on attenuation values (see Sec. 3.1). Moreover, the bone density also appears in the denominator of the mechanical stimulus of (7) discussed in the following paragraphs. To avoid a division by zero, we thus also assign ρ0=ρ~b for all x ∈ Ωs.

The rate of change in bone density ρ˙(x) is computed at each time-step as [27]
ρ˙(x)=SeffS(x)r˙(x)ρ^b
(4)
where S is the specific surface area in element e. Seff is an empirical constant that represents what fraction of the available surface area actually supports bone growth. S is calculated differently for the bone and scaffold domains and will be covered in subsequent sections. Likewise, Seff has different values for the two regions.
The bone density growth rate is calculated as [27]
r˙(x)={csn¯(x)(Ψ(x)Ψ*w)ifΨ(x)Ψ*>w0otherwise
(5)
where r˙ is the deposition rate of bone in units of length per unit time; n¯ is the normalized concentration of osteoblasts (with n¯(xΩb)=1); Ψ is the mechanical stimulus at element e; Ψ* is an empirical constant representing a reference stimulus level; w is an empirical constant representing the half-width of the dead zone around the reference level where neither bone growth nor resorption occur; and cs is an empirical constant.
The concentration in the scaffold, n¯(xΩs), is determined via a mass diffusion analysis by solving Fick’s second law
n¯t=μΔn¯
(6)
where μ is the diffusivity tensor and Δ is the Laplacian operator. We note that this equation has the same form as the heat equation and so is equivalent to a simplified unsteady heat conduction problem. The computation of the diffusivity tensor for the bone and scaffold regions is detailed in Secs. 2.2.1 and 2.2.2, respectively. The mechanical stimulus for N loadings is found via [11]
Ψ(x)=(ρ^bρt(x))2(i=1Nniσ¯i(x)m)1/m
(7)
where ni is the number of cycles of loading i the bone experiences per unit time, σ¯i is the effective stress for loading i, and m is an empirical weighting factor [11]. We note that the second factor on the right-hand side of this expression is similar to a p-norm that approximates the largest effective stress among all loading scenarios. Thus, this mechanical stimulus model assumes that the stimulus is driven by the largest effective stress (the true maximum is attained when m → ∞). The first factor, on the other hand, renders a larger stimulus for smaller values of the bone density ρt. Note that this factor also means the bone density can never be allowed to be exactly zero anywhere, else a division by zero will occur. We thus impose a minimum bone density everywhere, denoted as ρb~.
The effective stress is computed as [43]
σ¯(x)=2EAvg(x)W(x)
(8)
where W is the strain energy density. EAvg is the average of the material’s Young’s moduli. In the case of an isotropic material, it simplifies to the elastic modulus. The values for all empirical constants used herein are shown in Table 1.

2.2.1 Bone Region Update.

The elastic modulus and Poisson’s ratio in the bone region, x ∈ Ωb, are obtained via [41]
E(x)={2014(ρt(x))2.5,ρt(x)<1.2g/cm31763(ρt(x))3.2,ρt(x)1.2g/cm3
(9)
ν(x)={0.2,ρt(x)<1.2g/cm30.32,ρt(x)1.2g/cm3
(10)
where E is the elastic modulus in MPa, ν is the Poisson’s ratio, and ρt is the bone density. This model assumes that the elastic properties of bone are isotropic (albeit heterogeneous), as in Ref. [27]. The specific surface in the bone region is calculated using [44]
S(xΩb)=Sb(x)=32.3ϕ(x)93.9ϕ(x)2+134ϕ(x)3101ϕ(x)4+28.8ϕ(x)5
(11)
which has units of mm−1, with
ϕ(x)=1ρt(x)ρ^b
(12)
being the porosity of the bone. The effective fraction of specific surface area for the bone region is denoted as Seff(xΩb)=Seffb and its value is listed in Table 1. Note that in our implementation, ϕ never attains a value of 1 due to the lower bound on the bone density, namely, 0ϕϕmax=1ρ~b/ρ^b.
The diffusion properties of the bone are updated at every time-step for the mass diffusion analysis. The diffusivity of bone is estimated via [27]
D(x)=D^bk(x)k^b
(13)
where
k(x)=Cϕ3(x)S(x)2
(14)
is the permeability of bone as determined by the Kozeny equation, C is the Kozeny constant of bone, and D^b is the maximum diffusivity of bone (Table 1). The maximum permeability of bone, k^b, is computed using (14) and (11) by setting ϕ = ϕmax. Here too it is assumed that bone is isotropic, thus μ(x) = D(x)I in (6), where I is the identity two-tensor.

2.2.2 Scaffold Region Update.

Due to the arrangement of layers in the rectilinear scaffold considered here (shown in Fig. 1), where rods in alternating layers form a right angle, we expect the elasticity and diffusivity tensors of a homogenized material that behave in the bulk as the scaffold to be orthotropic. We detail in Sec. 2.3 the efficient computation of these homogenized properties. The homogenized moduli of the scaffold are used in (8) to compute the effective stress.

As we are using a homogenized material to model the scaffold, an assumption built in the update of (3) in the scaffold region is that when bone grows at a point x, it entirely fills in the scaffold at that point. That is, at time t, if ρt>ρ~b, the interstitial space at x does not contain void regions and is full of bone with density ρt. Note that this lower bound ρ~b is imposed for numerical reasons (cf. Sec. 2.2), since if the bone density is zero, then the mechanical stimulus (7) is undefined. In reality, when the scaffold is implanted, the bone density inside the scaffold is clearly zero. This assumption of a uniform bone distribution is reasonable for a small enough finite element size.

To account for the rotation θ of the scaffold about the axis perpendicular to the deposited layers in the bone adaptation analysis, the homogenized elasticity and diffusion tensors computed using the surrogates described in Sec. 2.3 are rotated according to
C^ijkl=RipRjqRkrRlsCpqrsH
(15)
μ^ij=RipRjqμpqH
(16)
respectively, where the summation convention is used. R is the rotation matrix given by
R=[cosθsinθ0sinθcosθ0001]
(17)
which assumes that the out-of-plane axis of the scaffold is the three-axis. This would be the z axis as shown in Fig. 3. Note that we are not rotating the unit cells relative to one another, but rather rotating the scaffold as a whole with all its unit cells. This is depicted in the fourth row of Fig. 3.
We estimate the specific surface area for the rectilinear scaffold region as
S(xΩs)=Ss(x)=ρ~bρt(x)Ss0+Sb(x)
(18)
where
Ss0=2πd(l2αd)V=π(ps+d(12α))(1α)l2
(19)
is an approximation of the specific surface of the empty scaffold, with V = 2al2 the volume of the unit cell and
ps=ld
(20)
the pore size. Equation (18) corresponds to an approximate interpolation between two extremes: the first corresponds to the absence of bone in the scaffold, in which case ϕ = ϕmax ≈ 1, Sb(x) ≈ 0 and ρt=ρ~b, and therefore SsSs0. The second case occurs when the scaffold is occupied by fully dense bone, ρt=ρ^b, in which case 0<ρ~b/ρ^b1 and consequently Ss(x) = Sb(x) ≈ 0. Therefore, as bone starts growing inside the scaffold, the available surface area for additional bone to grow starts decreasing.

The effective elastic and diffusive properties of the scaffold are independent of its length scale and dependent only on the dimensionless parameters d/l and α. However, as seen in (19), changes in the length scale (e.g., the pore size) affect the specific surface area S. Since the rate of change of bone density depends on S, cf. (4), the bone growth depends on the length scale. This is consistent with empirical evidence that certain pore size ranges favor bone growth in scaffolds [45]. When the bone has reached the maximum density ρ^b so that ϕ ≈ 1, the scaffold’s specific surface Ss is approximately zero and therefore the rate of change in bone density of (5) is approximately zero.

The effective fraction of the specific surface area for the scaffold region is given by Seff(xΩs)=Seffs, see Table 1. The diffusivity of the scaffold is directly obtained via homogenization, i.e., μ(x) = μH(x), see Sec. 2.3.

2.3 Scaffold Properties Surrogate Model.

To compute the expressions of Sec. 2.2.2, we construct a surrogate model of the effective elastic and diffusivity properties of the scaffold as functions of the geometric parameters. This allows for the efficient, off-line computation of the properties for any value of the scaffold design parameters during the bone growth simulations. The geometric variables of the scaffold used as inputs for the surrogate are d/l, α, and ps (in μm).

In addition to these variables, we also consider the effect of bone ingrowth. Thus, the bone density within the interstitial spaces of the scaffold (ρ) is also an input to the surrogate. This is an important feature of our model, since as bone grows into the scaffold in the bone growth simulation, it provides an increasing amount of mechanical support and reduces the diffusivity. The mechanical properties of the bone in the interstitial region of the unit cell are computed using the expressions of Sec. 2.2.1. In this aspect, our model improves on the approximation of effective properties used in Ref. [27], since in that work the effect of bone growing in the scaffold is only taken into account as a change in the local porosity and specific surface area of the scaffold, which indirectly assumes that the properties of bone growing inside the scaffold are the same as those of the scaffold material.

To build the scaffold properties surrogate, we perform numerical homogenization of the scaffold–bone periodic unit cell using finite element analysis (FEA) for a number of designs (i.e., combinations of values of d/l, α, and ρ). This method is commonly used to approximate the effective bulk properties of periodic lattice structures such as the scaffolds considered here [48]. It is worth noting that this method considers the material to be infinitely periodic and thus it ignores the actual shape of the entire scaffold, which would be manufactured to fit the shape of the bone defect. For instance, if the defect shape is cylindrical, as the one considered in the example of Sec. 4, then a cylinder-shaped scaffold such as the one shown in Fig. 1 would be manufactured. Nevertheless, since the size of the scaffold unit cell considered here is small relative to the size of the defect, we assume that the homogenized properties are a reasonable approximation of the effective properties at any point within the scaffold. This assumption is typically made in, for instance, homogenization-based topology optimization techniques for multi-scale design.

To efficiently generate the finite element model for these designs, we employ the geometry projection method [36,37]. This method projects the geometry of the unit cell onto a density field discretized on a fixed finite element grid. This allows us to circumvent the challenge and expense of generating a body-fitted mesh for every different design.

In the geometry projection method, geometric primitives described by high-level parameters are mapped onto a density field. This is then discretized onto a fixed mesh for analysis. Here, the rods in the scaffold unit cell are represented as cylinders, parameterized by the cylinder’s axis endpoint positions and its radius. The unit cell is meshed with a regular grid of hexahedral elements, see Fig. 5. The projected density ηqe corresponding to rod q for each element e in this mesh is computed as (see, for example, Ref. [38])
ηqe={1,dqerf12+14(dqerf)33dqe4rf,rf<dqe<rf0,dqerf
(21)
where dqe is the signed distance from the centroid of element e to the boundary of rod q, and rf is the radius of the smallest sphere that circumscribes the element. Equation (21) is a smooth approximation of the Heaviside function (cf. Ref. [38]). The reader is referred to Refs. [37,38] for expressions to compute the signed distance as a function of the rod geometric parameters.
Fig. 5
Geometry projection of unit cell. (a) Projected density ηe on unit cell mesh. (b) ηe = 0.5 iso-surface.
Fig. 5
Geometry projection of unit cell. (a) Projected density ηe on unit cell mesh. (b) ηe = 0.5 iso-surface.
Close modal
The projected density in regions where rods intersect is given by the maximum of the projected densities for the intersecting rods. We replace the non-differentiable maximum function with a smooth approximation. Here, we use the smooth maximum introduced in Ref. [39]:
ηe=q=0nbarH~(dqe,ε)ηqeq=0nrodH~(dqe,ε)+ξ
(22)
where ξ is a small positive number to avoid division by zero for elements outside the rods, nrod = 6 is the number of rods in the unit cell, and H~ is the smooth approximation of the Heaviside function given by
H~(x,ε)={1,x<ε12x2ε+12πsin(πxε),εxε0otherwise
(23)
The parameter ɛ in this expression corresponds to the width of the transition region (from 0 to 1) of the smooth Heaviside; here it is set equal to d/100. A reason to use this particular smooth maximum approximation is that ηe ≤ 1, which is not the case for other approximations (e.g., the p-norm employed in Ref. [49]), which may render values larger than 1 and subsequently produce negative effective diffusivity values.

The elements in the mesh are sized to keep them as close to cubic as possible while maintaining a minimum of 20 elements per side of the unit cell. An example of the geometry projection for a unit cell is shown in Fig. 5.

As in density-based topology optimization techniques, we use the projected density to define an element-uniform ersatz material with elasticity tensor given by
Ce=Cb(ρ)+ηe(CrCb(ρ))
(24)
where Cb and Cr are the elasticity tensors of bone (computed using the expressions in Sec. 2.2.1) and the material of the rods, respectively. We note that the values of the bone density ρ and therefore of Cb are the same for the entire unit cell. We also note that since a very small but positive bone density is specified in the scaffold in the initial conditions (for the reasons outlined in Sec. 2.2), a zero elasticity tensor is not possible and thus the analysis is well posed.
Similarly, the ersatz diffusivity tensor is given by
μi=μb(ρ)+ηe(μ0μb(ρ))
(25)
where μb is the diffusivity tensor of bone (computed using the expressions in Sec. 2.2.1) and μ0 = D0I is a non-zero but small diffusivity tensor to prevent an ill-posed analysis (since we assume that the material of the rods has zero diffusivity); the value of D0 is listed in Table 1.
To compute the effective elastic properties of the scaffold, we use asymptotic homogenization (cf., for example, Refs. [35,50]). For a given scaffold design, this requires six finite element analyses of the unit cell corresponding to six applied unit strains ɛ;0(ij) = eiej (with ei the unit vector along the ith coordinate) and with periodic boundary conditions imposed on all six faces. The components of the homogenized elasticity tensor are given by
CijklH=1VVCpqrs(εpq0(ij)εpq(ij))(εrs0(kl)εrs(kl))dv
(26)
where Cpqrs are the components of the ersatz elasticity tensor of (24), V is the unit cell volume, and ɛ;(ij) is the strain resulting from the analysis with applied ɛ;0(ij) [33].
Similarly, to perform the homogenization of the diffusivity tensor, three unit concentration gradients c0(i)=ei are imposed, and the components of the homogenized tensor are given by
μijH=1VVμlm(c,l0(i)c,l(i))(c,m0(j)c,m(j))dV
(27)
where μlm are the components of the ersatz diffusivity tensor of (25), and c(i) is the concentration gradient resulting from the analysis with applied c0(i) [33].

The numerical implementation of the homogenization analysis was done in matlab version R2019b, which is based on the implementations of Refs. [33,34]. The rods are assumed to be made of a hydroxyapapatite/tricalcium phosphate (HA/TCP) mixture, with an elastic modulus of 10 GPa and Poisson’s ratio of 0.25 [27].

To build a surrogate model, a parametric sweep was conducted for bone density ρ{0.05,0.06,,1.85,ρ^b} g/cm3, rod overlap α ∈ {0.05.0.1, …, 0.45}, and rod diameter to separation ratio d/l ∈ {0.2, 0.3, …, 0.8}, for a total of 1260 combinations. The d/l and α ranges are derived from known limitations in manufacturing. The pore size limits are based on the range of pore sizes over which bone has been found to grow best in a scaffold [51] Accordingly, the surrogate model should only be considered valid in this range. The surrogate model employed for the scaffold–bone effective properties is a piecewise-multilinear interpolation [52]. This model was chosen due to its simplicity and to avoid the oscillations other models exhibited between reference points. Plots of some of the homogenized properties are shown in Fig. 6.

Fig. 6
Plots of some of the homogenized scaffold properties. The plots on the left correspond to the properties of an empty scaffold; the ones on the right correspond to the properties of a scaffold filled with fully dense bone. In-plane and out-of-plane are as defined in Fig. 2. (a) In-plane elastic modulus. (b) Out-of-plane elastic modulus. (c) In-plane diffusivity. (d) Unit cell axis reference.
Fig. 6
Plots of some of the homogenized scaffold properties. The plots on the left correspond to the properties of an empty scaffold; the ones on the right correspond to the properties of a scaffold filled with fully dense bone. In-plane and out-of-plane are as defined in Fig. 2. (a) In-plane elastic modulus. (b) Out-of-plane elastic modulus. (c) In-plane diffusivity. (d) Unit cell axis reference.
Close modal

2.4 Bone Growth Surrogate Model.

The purpose of the bone growth surrogate model is to efficiently approximate the quantities of interest in the scaffold design (e.g., the stiffness at the time of implantation and the bone growth at observation time) as a function of the scaffold design parameters. The mechanobiological simulation is computationally expensive, with run times typically measured in hours. Although analytical sensitivities could be employed for the optimization, their derivation and computational implementation are quite involved, particularly since this is a time marching problem. Furthermore, this would require smoothing of some of the non-differentiable functions in the bone growth model, for instance, the bone density growth rate of (5) and the bone mechanical properties of (9) and (10). Since there are only four design parameters (namely, α, d/l, θ, and ps), we therefore choose to employ a surrogate model here and defer the computation of analytical sensitivities for future work. Another advantage of using the surrogate is that we can choose a priori the number of computational bone adaptation simulations to run, and these simulations can be performed in parallel. The optimization is subsequently performed using the surrogate model. As customary, a full-fledged bone adaptation simulation of the optimal scaffold found with the surrogate is then performed.

To create the surrogate model, bone adaptation simulations were performed for a Latin hypercube sampling (LHS) of the design space. The ranges α ∈ [0.05, 0.45], d/l ∈ [0.2, 0.8], ps ∈ [100, 500] μm, and θ ∈ [0, π/2] rad were used. Note that, because of the rods in alternating layers of the rectilinear scaffold are perpendicular, it is only necessary to consider an orientation θ within one quadrant.

To select the surrogate model to fit the data, k-fold cross validation was employed (cf. Ref. [53]). This is performed by first randomizing the order of the simulation results, then subdividing the results into subgroups (or folds) of k members each. The surrogate model is fitted to all of the data minus one of the subgroups. The inputs of the subgroup left out are then plugged into the fitted surrogate model, and the predicted results are compared with the known output values from the subgroup. This process is repeated for every subgroup [54]. The suitability of a surrogate model is determined based on how well the surrogate performs across all subgroups [53]. In our case, we chose the model that exhibits the highest minimum R2 value across all folds. As detailed in Sec. 4.2, the surrogate model selected based on this procedure for the example presented in this paper is Kriging.

We construct surrogates of three quantities from the bone adaptation simulation. The first is the strain energy of the scaffold at time of implantation (t0), given by
U0=12Ωs[σ(x):ε(x)]|t=0dveΩsse|t=0
(28)
where se is the element strain energy. This quantity is equivalent to half of the structural compliance of the scaffold and is an inverse measure of stiffness. It is important because it relates to the ability of the scaffold to sustain loads upon implantation.
The other quantities approximated via surrogates relate to bone growth. The first corresponds to mass bone growth fraction at the final simulated time (tf), given by
mf:=Ωsρt(x)|t=tfdvΩsρ^bdveΩsρe|t=tfveeρ^bve
(29)
where ρe is the bone density of element e and ve is the element volume. A value of mf = 1 is obtained if the scaffold is full of bone at its maximum density (ρ=ρ^b at the final simulation time). Henceforth, we will refer to t0 as the “implantation time” and tf as the “final time.” This fraction is commonly used to quantify bone growth in the literature (see, for example, Ref. [27]).
One possible way to define what complete osteointegration means is to have a mass bone growth fraction of 1 everywhere in the scaffold. In this situation, fully dense bone occupies the entire scaffold. However, a more realistic definition may be to require that bone occupies the entire scaffold, but not necessarily at the maximum density. The latter definition is more consistent with the fact that the bone itself is not fully dense everywhere. The majority of reported in vivo studies have shown insufficient bone growth, with growth only in the scaffold periphery, which has been postulated to be due to insufficient vascularization [55,56]. This lack of vascularization leads to a shortage of oxygen and glucose needed for viable tissue and causes tissue necrosis. To better capture the ability of the bone to occupy all of the interstitial space in the scaffold regardless of density, we define a volume bone growth fraction at the final simulated time as
vf:=ΩsH(ρt(x)ρthr)|t=tfdvΩsdveΩsH(ρeρthr)|t=tfveeve
(30)
where H is the Heaviside function H(x)={0ifx<0,1otherwise}, whose role is to ensure that only elements with a bone density equal to or greater than a specified threshold bone density ρthr are counted in the sum. That is, vf is a volume fraction of the interstitial space in the scaffold that has some bone growth at the final simulation time. How exactly ρthr is selected is covered in more detail in Sec. 3.1. In Sec. 4, we report these bone growth fractions as percentages, Gv = 100 vf and Gm = 100 mf.

As detailed in Sec. 4.3, we employ these surrogates to optimize the scaffold design. The goal of the optimization is ensuring to maximize the bone growth at the final simulation time while ensuring that the scaffold possesses enough stiffness at implantation time.

3 Implementation

To demonstrate our methodology, we consider a cylindrical scaffold implanted in an adult rabbit femur. The scaffold has a 6 mm diameter and it goes through the bone perpendicular to the frontal plane (approximately 10 mm in length).

3.1 Methods.

Micro-CT images of a rabbit femur were acquired from the harvested femur of a male white New Zealand rabbit weighing approximately 3.5 kg. Prior to euthanization, the animal was housed in an individual cage with ad libitum access to food and water. Care and use of the animal are adhered to University of Connecticut’s Institutional Animal Care and Use Committee guidelines. The left femur of the rabbit was dissected and stored fresh frozen until the time of the scan. For micro-CT imaging, the femur was bisected with a rotary tool at the mid-diaphyseal region. The proximal portion was discarded, and only the distal portion is modeled in the simulation. The specimen was scanned using a μCT50 (Scanco Medical, Switzerland) after positioning it vertically in the cylindrical container. The scanning protocol included 1000 projections, an energy of 55 V, an intensity of 145 μA, and a voxel size of 34.4 μm. The geometry was extracted using 3d slicer 4.10.2. The hole where the scaffold is implanted was created in the imported geometry and not in the physical sample. Computer-aided design (CAD) models of the femur and scaffold regions created in 3d slicer are shown in Figs. 7 and 8, respectively. The facets shown in these figures do not correspond to the finite element mesh used for the analysis but to the tessellated boundary representation created from the micro-CT scan.

Fig. 7
Femur geometry with defect, no scaffold. The femur geometry is derived from the micro-CT scan. (a) XY view and (b) YZ view.
Fig. 7
Femur geometry with defect, no scaffold. The femur geometry is derived from the micro-CT scan. (a) XY view and (b) YZ view.
Close modal
Fig. 8
Scaffold region/bone defect. This is the region the scaffold is expected to occupy and into which we hope to induce bone growth.
Fig. 8
Scaffold region/bone defect. This is the region the scaffold is expected to occupy and into which we hope to induce bone growth.
Close modal

It is important to note that the exact geometry of the scaffold itself (i.e., the rods) is not modeled in the mechanobiological simulation. Instead, this region is assigned a homogenized material with properties equivalent to those of a scaffold composed of a specified unit cell design filled entirely with bone at a given density. These properties are derived from the scaffold surrogate (cf. Sec. 2.3), and they are updated on an element-by-element basis at each time-step of the transient mechanobiological simulation. Thus, the rectilinear scaffold structure is not visible in Fig. 8.

The initial bone density distribution for the femur was derived from the micro-CT scan data. Figure 9(a) shows a slice of the femur scan. Micro-CT scans are effectively 3D maps of the x-ray attenuation of the sample, typically recorded in Hounsfield units [57]. The attenuation of bone has been found to be directly proportional to its density [58,59]. Thus, we may calculate the bone density from the attenuation values in the micro-CT scan as ρvoxel = voxel + b, where m and b are empirical constants specific to a given micro-CT scanner and combination of scan settings, and ζvoxel is the attenuation value. Ordinarily, m and b are found using calibration phantoms. However, the micro-CT scan was not furnished with any calibration information, so we approximated m and b as m=ρ^b/max(ζ) and b = 0, where ζ is the vector of all attenuation values in the scan. These values were then mapped to the mesh of the femur by assigning to each element the density of the micro-CT scan voxel closest to the element’s centroid. The density values mapped in this way onto the meshed femur corresponding to the slice of Fig. 9(a) are shown in Fig. 9(b).

Fig. 9
Micro-CT scan slice and corresponding bone densities mapped to the FEA mesh. (a) CT scan and (b) FEA femur model.
Fig. 9
Micro-CT scan slice and corresponding bone densities mapped to the FEA mesh. (a) CT scan and (b) FEA femur model.
Close modal

We also use the initial density field obtained from the micro-CT scan to determine a reasonable value for the bone density threshold ρthr of (30), below which we consider element to be primarily unmineralized soft tissue that should not be counted towards bone growth. This visual selection of a threshold was done because there is no widely agreed upon density cutoff value for what should be considered bone in a micro-CT scan [60]. Indeed, the cutoff values in some studies are determined by visual inspection of the micro-CT scan as we did [61,62]. A variety of different values have been used in micro-CT imaging studies, typically given either as a percentage of the maximum attenuation in the scan or in mg/ml of hydroxyapatite (HPA) [60]. The latter can be interpreted as the attenuation of a calibration phantom that contains concentration of HPA, a common type of calibration phantom for bone. Example thresholds in literature vary from 320 mg/ml HPA to 960 mg/ml HPA, and from 11% to 54%. This corresponds to a bone density range of 0.2 g/cm3 to >1.6 g/cm3 [60,63]. Thus, we chose ρthr to be approximately equal to the trabecular bone density in the region where the scaffold is implanted, which is roughly 0.75 g/cm3. This is within the range of cutoff values reported for studies targeting trabecular bone, equivalent to a threshold of 39%.

3.2 Implementation of Bone Adaptation Analysis.

Both diffusion and mechanical analyses were performed using abaqus [64]. We employed user subroutines to ensure that the time stepping for a single bone adaptation simulation was done in a single run, as opposed to invoking a new abaqus analysis at each time-step to avoid the overhead and make the simulation more efficient. The pseudocode for the simulation procedure is presented in Fig. 10.

Fig. 10
Pseudocode for the mechanobiological analysis loop
Fig. 10
Pseudocode for the mechanobiological analysis loop
Close modal

To perform the mechanical and diffusion analyses simultaneously, we employed the combined heat transfer-static analysis in abaqus. While the solver does support mass diffusion analysis, the required element types are not compatible with those of a static analysis, so they could not both be accomplished in a single run. Moreover, it is not possible to change diffusion properties while the simulation is running, which is necessary since the effective elastic and diffusive properties of the bone–scaffold system change over time. As noted in Sec. 2.2, with the specific heat and material density set to unity, a heat transfer analysis is identical to a diffusion analysis, and the user subroutine UMATH allows to set the properties in between time-steps. The thermal conductivity of the material was thus used as the analog for diffusivity and temperature for concentration. The UEXTERNALDB subroutine was used to pause the analysis between load steps to calculate the new material properties, and the UMAT and UMATH subroutines were used to update the material properties in the analysis.

The ORIENT subroutine was used to compute the rotation of the scaffold elasticity and conductivity tensors arising from the scaffold orientation angle θ corresponding to (15) and (16), respectively. The calculation of the effective properties was performed using a python script.

A single loading case corresponding to hopping was considered, with an estimated force magnitude of 184 N [65], and a load frequency of 6000 cycles per 6 h applied to the femur for the duration of the simulation. The simulated time period starts at the completion of the implantation procedure or implantation time (t0 = 0 days) and ends at 28 days later (tf = 28 days) at the final time [27]. The simulation thus comprises 113 time-steps (including the initial state). The force was split between the two condyles and applied as a surface traction spread over several nodes on each condyle. The tractions are parallel to the sagittal plane, and they form a 55 deg angle with respect to a line perpendicular to the femur’s transverse plane. Zero-displacement boundary conditions were applied at the cut surface of the femur. The boundary conditions for the mechanical analysis are shown in Fig. 11.

Fig. 11
Femur mechanical boundary conditions. Again, as in Figs. 7 and 8, no mesh is shown. The arrows indicate the direction and region of application of the knee joint load to the femur. The triangular markers at the cut plane of the femur indicate an encastre boundary condition across its surface. (a) XY view and (b) YZ view.
Fig. 11
Femur mechanical boundary conditions. Again, as in Figs. 7 and 8, no mesh is shown. The arrows indicate the direction and region of application of the knee joint load to the femur. The triangular markers at the cut plane of the femur indicate an encastre boundary condition across its surface. (a) XY view and (b) YZ view.
Close modal

For the diffusion analysis, the bone of the femur surrounding the scaffold is treated as a source of osteoblasts, with a fixed normalized concentration of 1. The scaffold is assumed to have a zero initial osteoblast concentration at implantation [27]. The coefficient of thermal expansion is set to zero for both bone and scaffold to ensure the temperature, which is the proxy quantity for concentration, does not affect the mechanical analysis.

The femur and scaffold are meshed with 225,022 linear tetrahedral elements for coupled temperature–displacement analysis. The meshed model is shown in Fig. 12.

Fig. 12
The meshed femur and scaffold assembly, showing elements in the bone region Ωb and in the scaffold region Ωs (cylindrical region towards distal end of femur).
Fig. 12
The meshed femur and scaffold assembly, showing elements in the bone region Ωb and in the scaffold region Ωs (cylindrical region towards distal end of femur).
Close modal

4 Results

In the following, we describe the numerical design of experiments employed to build the surrogate models of bone growth and strain energy, and the optimization of the scaffold based on these surrogates.

4.1 Latin Hypercube Sampling.

We performed a full bone adaptation simulation for the 81 designs of the LHS. Scatter plots of volume and mass bone growth and strain energy for these runs are shown in Fig. 13. These plots show the combined effect of two variables on the responses of interest. Where multiple markers may overlap, the one attaining the largest bone growth value is plotted in the foreground.

Fig. 13
Scatter plots of bone growth and strain energy for the scaffolds in the LHS. Each marker is colored based on its (a) volume or (b) mass growth percentage at t = tf. Marker size indicates strain energy of the scaffold at implantation. (a) Volume growth and (b) mass growth.
Fig. 13
Scatter plots of bone growth and strain energy for the scaffolds in the LHS. Each marker is colored based on its (a) volume or (b) mass growth percentage at t = tf. Marker size indicates strain energy of the scaffold at implantation. (a) Volume growth and (b) mass growth.
Close modal

Several important trends appear in the results. From the perspective of the scaffold performance, we aim to have maximum bone growth (i.e., red markers) at the final time for a given strain energy at implantation time. The red markers are the largest across all the d/l plots of Fig. 13, which indicates that there is a trade-off between either of the bone growth measures and stiffness. Figure 13 also shows that d/l has the largest effect on bone growth. Values of d/l approaching the lower bound of 0.2 show the greatest bone growth as measured by either volume or mass. Conversely, values of d/l approaching the upper bound of 0.8 render the highest stiffness.

Another interesting observation is that while pore size has some effect on bone growth, it is not as significant of a driver as d/l. Indeed, varying the pore size along the bounds ps ∈ [100, 500] μm for the stiffest scaffold with d/l = 0.8 does not render as significant an improvement in bone growth as compared to simply decreasing d/l, which is an indication that stiffness is the major driver of bone growth. Pore size still has some effect; this is expected because the specific surface area of (4) does depend on the actual dimensions of the scaffold, with smaller pores yielding a larger specific surface area. The bone thus has a larger surface to grow on and less volume to fill. However, decreasing pore size necessarily reduces the diffusivity of the scaffold, making the growth front harder to reach by the biological actors involved in bone growth [21,27,51]. We do not see much of an effect here, however, likely because the range of pore sizes was deliberately chosen for good bone growth as mentioned in Sec. 2.3.

It is notable that the simulation results indicate there is some freedom to select d or l. As the d/l ratio is unitless, the designer may choose either (though not both) to suit the application without compromising bone growth. This observation has important manufacturing and mechanical repercussions. For instance, a rod diameter that is more amenable to the manufacturing process (e.g., one that renders less variability in dimensions and produces less manufacturing defects) can be selected largely independently of the impact on bone growth. Moreover, a larger d requires fewer rods in a scaffold of the same size as compared to one with smaller d for the same d/l, which would decrease the deposition time. Another consideration is that rods with smaller d may exhibit a larger compressive strength due to size effects in the case of brittle materials (such as HA/TCP). Therefore, a scaffold with smaller rods (but the same overall volume fraction) may be stronger in compression. Moreover, if there are more rods in the scaffold, there is more redundancy should a rod break, and thus the scaffold would be better from the point of view of fail-safe design.

Another important observation from the simulations is that large values of ps may hinder growth. Since all the other defining geometric variables are dimensionless in our parameterization, pore size establishes the length scale of the unit cell. Per the square-cube law, the surface area increases quadratically and volume cubically as the scale of an object increases. Therefore, the specific surface area of the empty scaffold Ss0 in (19) is inversely proportional to pore size (assuming all other dimensions are held constant). Thus, if ps is large then the specific surface area S for the scaffold region of (18) will be small (since Sb ≈ 0 in the empty scaffold). Consequently, the rate of change in bone density of (4) will be reduced.

Finally, we note that there seems to be little difference between the interaction plots for the mass and the volume bone growth fractions in Fig. 13. This is an indication that, at least for the problem considered in this work, larger bone mass equates to larger occupied volume. We note, however, that this may not be the case in other situations, where fully dense bone may quickly form around the periphery of the scaffold and subsequently stop growing towards the center of the scaffold. In that case, the mass growth fraction may be larger than the volume growth fraction, since the latter can be considered an inverse measure of void space in the scaffold with no or limited bone growth. We should note, however, that there is a more noticeable difference between the two measures when considering the effect of individual parameters, as discussed in the following section.

We postpone our observations on the effect of the overlap ratio α and the orientation θ on bone growth to the following section, where their influence can more clearly be observed.

4.2 Surrogate Model.

A surrogate model was fitted to the runs of the LHS for all performance metrics (U0, Gm, and Gv) using the Surrogate Modeling Toolbox (SMT) for python [66]. Using k-fold cross validation, a Kriging model was chosen for U0 and Gm, and a Kriging partial least-square model was used for Gv. The default parameters were used, with constant regression functions, squared exponential correlation functions, and a hyperparameter value of 0.01. These models rendered worst-case R2 values of 0.998, 0.997, and 0.949, respectively.

In Fig. 14, we can see a two-parameter contour plot of a 20-level, 4-factor full-factorial sampling of the surrogates. The color in each plot of the figure with parameters x1 and x2 in the axes corresponds to the maximum bone growth fraction that can be obtained with x1 and x2 fixed and the remaining two parameters being design variables.

Fig. 14
Filled contour plots of scaffold bone growth as predicted by the surrogate model. (a) Volume growth and (b) mass growth.
Fig. 14
Filled contour plots of scaffold bone growth as predicted by the surrogate model. (a) Volume growth and (b) mass growth.
Close modal

Figure 14 confirms the previously made observations about the effect of d/l and ps on bone growth (which is expected since the surrogate is built on the LHS data). However, two additional observations become more evident in this figure.

The first is that the orientation θ has almost no effect on bone growth. This occurs because the largest bone growth occurs in the scaffolds with lower d/l, and for these scaffolds it can be readily shown via the surrogates that the magnitude and variation of Eavg in (8) with respect to θ are relatively small for most of the angle range, with higher values attained around θ = 0 (which is equal to those for θ = π/2). While the insensitivity of bone growth to orientation may hold for the particular defect geometry, loading and boundary conditions considered in our example, it most likely cannot be generalized to other cases. It is possible that higher sensitivity could be observed in other situations—e.g., a defect in a flat bone situated in a region loaded primarily in tension.

The second observation is the effect of the overlap fraction α on bone growth. Although its effect is clearly not as important as that of d/l, it can be observed from Fig. 14 that the smallest α, the largest the bone growth. This makes sense from a stiffness point of view, since decreasing α decreases the stiffness of the scaffold, and as the strong influence of d/l on growth indicates, a more compliant scaffold fosters higher growth.

In the surrogate for volume growth of Fig. 14, there is some jaggedness in the θ versus d/l plot as d/l → 0.8. It is unclear what causes this phenomenon, but it is possible it may be the result of oscillations in the surrogate model between training points. Additionally, the ps versus d/l plot displays a region of slightly elevated volume growth around d/l = 0.65 and ps = 450 μm. This seems to indicate a beneficial interaction between the parameters which is not readily apparent in Fig. 13(a).

To further examine how each individual parameter affects bone growth, we plot the largest bone growth value that can be obtained for the parameter, as shown in Fig. 15. The curves in Fig. 15 can thus be interpreted as growth limit curves, i.e., the upper bound each variable places on bone growth. Each point in these curves represents the maximum achievable bone growth if the scaffold is optimized with the parameter corresponding to the curve held fixed and all other parameters as design variables. Any point (x, y) in each curve is the result of an optimization run performed using the surrogates, where y is the maximal bone growth obtained by fixing the value of the parameter of the curve to x and varying all other parameters. For example, each point in the curve corresponding to d/l is the result of a maximization of growth with the value of d/l fixed and with α, ps, and θ as design variables. All parameters are normalized as z^=(zzmin)/(zmaxzmin), where [zmin, zmax] is the range of the parameter (cf. Sec. 2.4). Each optimization was run with a basin-hopping algorithm with a Newton conjugate gradient method as an local optimizer. Both implementations were from the SciPy library [67].

Fig. 15
Maximum-growth or growth limit curves for each parameter. Parameters are normalized by their ranges considered in the LHS sampling. (a) Volume growth and (b) mass growth.
Fig. 15
Maximum-growth or growth limit curves for each parameter. Parameters are normalized by their ranges considered in the LHS sampling. (a) Volume growth and (b) mass growth.
Close modal

As previously observed, d/l is by far the most influential design variable, while α and ps have a more moderate effect. As already concluded from the previous discussion, θ seems to have the least influence, especially for mass growth. Unlike the interaction plots of Figs. 13 and 14, the curves of Fig. 15 reveal differences between the mass and volume growth fractions. First, d/l and α have a more pronounced effect on the volume growth fraction than on the mass growth fraction. These curves also show that while the maximum possible volume and mass growth fractions are similar, the former varies over a greater range. This greater range of growth values indicates there is more volume growth to be gained by optimizing the scaffold design. Finally, the d/l limit curve for volume growth is not monotonic, which is likely due to the arbitrary choice of the bone density threshold ρthr of (30). The lack of monotonicity of volume growth with respect to the design variables justifies the use of a non-gradient-based method for the exploration stage of the optimization (a gradient-based method is subsequently used for exploitation). Interestingly, the d/l curve is non-monotonic for the mass growth as well. This is likely due to an increased specific surface area at high d/l due to (20) partially compensating for their high stiffness.

As noted in Sec. 2.4, in addition to trying to maximize bone growth, we also must ensure that the scaffold has sufficient stiffness to sustain the loads at the time of implantation. Figure 16 shows the Pareto front between the bone growth measures and the strain energy. The Pareto fronts show a clear trade-off between strain energy and growth. In general, strain energy is proportional to growth. This makes sense since bone growth is stimulated by mechanical strain. Comparing the LHS simulation runs to the Pareto front, we observe that even at fixed strain energy, there is considerable latitude to improve bone growth as measured by either volume or mass bone growth. An agglomeration of data points is present along the vertical axis of the volume growth measure Pareto front (Fig. 16(a)). This is due to designs where the bone growth never exceeds the stipulated threshold bone density ρthr.

Fig. 16
Pareto front of bone growth versus strain energy. The front (circular markers) was generated using the surrogate model, with runs in the LHS sample (cross markers) overlaid for reference. (a) Volume growth and (b) mass growth.
Fig. 16
Pareto front of bone growth versus strain energy. The front (circular markers) was generated using the surrogate model, with runs in the LHS sample (cross markers) overlaid for reference. (a) Volume growth and (b) mass growth.
Close modal

4.3 Optimization.

We now seek to demonstrate that the surrogate models can be used to improve the performance of the scaffold. Unfortunately, there are no in vivo rabbit studies that report bone growth in rectilinear scaffolds. Such research would provide a comparison with experimental data. There are, on the other hand, in vivo studies of this type of scaffolds implanted in pig mandibles (cf. Refs. [32,68,69]). Another study examined a scaffold of this type in sheep [30]. However, several of the model parameters and the loadings are not presently available for either sheep or pigs, and therefore we cannot apply our model to these cases. A comparison to experimental data is therefore postponed to future work.

As we lack experimental data for a proper comparison, we instead demonstrate the efficacy of the proposed design procedure by comparing the optimized scaffold to a reference scaffold, which we choose to have the same dimensions as the scaffolds used in the in vivo study of Lan Levengood et al. [32]. That is, we assume that the bone adaptation model is correct (an assumption that is reasonable in light of the comparison of this model to experimental data made in Ref. [27]), and we only examine the extent of the improvement that can be attained with the proposed design methodology. The in vivo study of Lan Levengood et al. [32] considers scaffolds implanted in pig mandibles; since the magnitudes of the loads experienced by the scaffold in the pig mandible are significantly larger than those of the scaffold in the rabbit femur, it is reasonable to assume that the reference scaffold has sufficient stiffness in the latter setting. We also note that the parameter values of the reference scaffold (listed in Table 2) are not too far from the midpoint values of the ranges considered in the LHS for each parameter.

Table 2

Comparison of the design and performance of the reference and optimized scaffolds

Scaffoldαd/lps (μm)θ (rad)U0 (N mm)Gv (%)Gm (%)
Reference0.180.52335902.348.222.8
Optimized0.450.34410002.291836.3
Scaffoldαd/lps (μm)θ (rad)U0 (N mm)Gv (%)Gm (%)
Reference0.180.52335902.348.222.8
Optimized0.450.34410002.291836.3

We use the surrogate models to find the scaffold design that maximizes the volume growth fraction subject to the constraint that the stiffness at the time of implantation is no smaller than that of the reference scaffold. The parameters derived from the surrogate model optimization as well as those for the reference scaffold were then evaluated using the full mechanobiological simulation. The performance of the optimized and reference scaffolds is listed in Table 2. As before, the surrogate optimization is performed using a basin-hopping algorithm, this time with a constrained trust region method as the local optimizer. The change in local optimizer is required by the addition of a stiffness constraint. Both are part of the SciPy library [67]. The performance predicted by the surrogate model shows a reasonable agreement with the values found via the full simulation. The estimated mass and volume bone growth were 15.3% and 34.6%, respectively, compared to the simulation derived 18% and 36.6%. The bone density of both designs at the final time tf is shown in Fig. 17.

Fig. 17
Bone density distribution of the reference (left) and optimal (right) scaffold designs at the end of the simulated time. Rows corresponding to different views. Only bone density is plotted, i.e., the scaffold density is not taken into account in the plots.
Fig. 17
Bone density distribution of the reference (left) and optimal (right) scaffold designs at the end of the simulated time. Rows corresponding to different views. Only bone density is plotted, i.e., the scaffold density is not taken into account in the plots.
Close modal

The optimized design is slightly stiffer than the reference design with more than 50% increase in mass growth fraction and nearly 120% increase in volume growth fraction. That the two scaffolds have similar stiffnesses but the optimized scaffold exhibits greater bone growth highlights an important result: optimizing the scaffold with regard to stiffness as a bulk material, which is the approach taken in the vast majority of scaffold design strategies published to date, does not necessarily guarantee maximum bone growth. The optimized scaffold also has a larger terminal rate of bone growth (Gv˙(t=tf) and Gm˙(t=tf)), as shown in Fig. 18. Neither scaffold has reached a steady state at the end of the simulation period, therefore it is possible that the optimized scaffold will reach even greater bone growth after the 28-day period considered in our simulation.

Fig. 18
Comparison of the growth histories for the optimized and reference scaffolds
Fig. 18
Comparison of the growth histories for the optimized and reference scaffolds
Close modal

5 Conclusions

Using a full mechanobiological simulation of bone growth, we optimized a rectilinear scaffold to maximize osteointegration. This growth-optimized scaffold exhibits greater bone growth when compared to a reference scaffold design that has been used in published in vivo studies (albeit in a different animal). Although this work presents a single case study, it is evident from the performance of the optimized design that incorporating bone adaptation directly into the optimization and designing the scaffold as a structure (i.e., taking into account the shape of the defect and the loading) can lead to significant improvements in bone growth. These findings require in vivo validation, which is the subject of ongoing work.

The proposed design methodology can be readily applied to other scaffold configurations and manufacturing processes. It has no theoretical restriction on the number of design variables, so it may be applied to any scaffold which can be geometrically parameterized and homogenized. The latter condition requires periodicity of the scaffold structure, which is a limitation of the method. However, there is no restriction on the shape and size of the scaffold that may be considered. The method also has a computational constraint; the number of simulations required to fit the bone growth surrogate grows exponentially with the number of design parameters, and so in practice the method is restricted by the available computational resources.

Several mechanobiological aspects of bone–scaffold design were not considered in this study and will be considered in future research. These include considerations such as stress and fatigue life of the bone–scaffold system and bone resorption. Moreover, modern manufacturing techniques afford far more flexibility in scaffold architecture than the rectilinear configuration examined here. Incorporating these considerations will expand the design space and thus potentially render better osteointegration, and they will increase the clinical viability of these scaffolds.

Acknowledgment

Support from the National Science Foundation, award CMMI-1727591, to conduct this work is gratefully acknowledged. We are also grateful to Dr. Yusuf Khan at the University of Connecticut Health Center for his help in obtaining the tissue for this study and preparing the tissue for scanning. We would also like to thank Prof. Mariana Kersh at the University of Illinois at Urbana-Champaign for useful discussions on this subject. Finally, we would like to thank Samuel N. Ferrer-Fuenmayor for performing some preliminary coding for the bone adaptation model.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtained from the corresponding author upon reasonable request.

References

1.
Finkemeier
,
C. G.
,
2002
, “
Bone-Grafting and Bone-Graft Substitutes
,”
J. Bone Joint Surg.
,
84
(
3
), pp.
454
464
.
2.
Dimitriou
,
R.
,
Jones
,
E.
,
McGonagle
,
D.
, and
Giannoudis
,
P. V.
,
2011
, “
Bone Regeneration: Current Concepts and Future Directions
,”
BMC Med.
,
9
(
66
), pp.
1
10
.
3.
Fernandez de Grado
,
G.
,
Keller
,
L.
,
Idoux-Gillet
,
Y.
,
Wagner
,
Q.
,
Musset
,
A. M.
,
Benkirane-Jessel
,
N.
,
Bornert
,
F.
, and
Offner
,
D.
,
2018
, “
Bone Substitutes: A Review of Their Characteristics, Clinical Use, and Perspectives for Large Bone Defects Management
,”
J. Tissue Eng.
,
9
.
4.
Khan
,
S. N.
,
Cammisa
,
F. P.
,
Sandhu
,
H. S.
,
Diwan
,
A. D.
,
Girardi
,
F. P.
, and
Lane
,
J. M.
,
2005
, “
The Biology of Bone Grafting.
,”
J. Am. Acad. Orthop. Surg.
,
13
(
1
), pp.
77
86
.
5.
Campana
,
V.
,
Milano
,
G.
,
Pagano
,
E.
,
Barba
,
M.
,
Cicione
,
C.
,
Salonna
,
G.
,
Lattanzi
,
W.
, and
Logroscino
,
G.
,
2014
, “
Bone Substitutes in Orthopaedic Surgery: From Basic Science to Clinical Practice
,”
J. Mater. Sci.: Mater. Med.
,
25
(
10
), pp.
2445
2461
.
6.
Brydone
,
A. S.
,
Meek
,
D.
, and
MacLaine
,
S.
,
2010
, “
Bone Grafting, Orthopaedic Biomaterials, and the Clinical Need for Bone Engineering
,”
Proc. Inst. Mech. Eng.
,
224
(
12
), pp.
1329
1343
.
7.
Wang
,
W.
, and
Yeung
,
K. W.
,
2017
, “
Bone Grafts and Biomaterials Substitutes for Bone Defect Repair: A Review
,”
Bioact. Mater.
,
2
(
4
), pp.
224
247
.
8.
Lauthe
,
O.
,
Soubeyrand
,
M.
,
Babinet
,
A.
,
Dumaine
,
V.
,
Anract
,
P.
, and
Biau
,
D. J.
,
2018
, “
The Indications and Donor-Site Morbidity of Tibial Cortical Strut Autografts in the Management of Defects in Long Bones
,”
Bone Joint J.
,
100-B
(
5
), pp.
667
674
.
9.
Metz
,
C.
,
Duda
,
G. N.
, and
Checa
,
S.
,
2019
, “
Towards Multi-Dynamic Mechano-Biological Optimization of 3D-Printed Scaffolds to Foster Bone Regeneration.
,”
Acta Biomater.
,
101
, pp.
117
127
.
10.
Wang
,
Y.
,
Arabnejad
,
S.
,
Tanzer
,
M.
, and
Pasini
,
D.
,
2018
, “
Hip Implant Design With Three-Dimensional Porous Architecture of Optimized Graded Density
,”
ASME J. Mech. Des.
,
140
(
11
), p.
111406
.
11.
Carter
,
D. R.
, and
Beaupré
,
G. S.
,
2001
,
Skeletal Function and Form: Mechanobiology of Skeletal Development, Aging, and Regeneration
,
Cambridge University Press
,
New York
.
12.
Adachi
,
T.
,
Osako
,
Y.
,
Tanaka
,
M.
,
Hojo
,
M.
, and
Hollister
,
S. J.
,
2006
, “
Framework for Optimal Design of Porous Scaffold Microstructure by Computational Simulation of Bone Regeneration
,”
Biomaterials
,
27
(
21
), pp.
3964
3972
.
13.
Christen
,
P.
,
Ito
,
K.
,
Ellouz
,
R.
,
Boutroy
,
S.
,
Sornay-Rendu
,
E.
,
Chapurlat
,
R. D.
, and
Van Rietbergen
,
B.
,
2014
, “
Bone Remodelling in Humans is Load-Driven But Not Lazy
,”
Nat. Commun.
,
5
(
4855
), pp.
1
5
.
14.
Egan
,
P. F.
,
Bauer
,
I.
,
Shea
,
K.
, and
Ferguson
,
S. J.
,
2019
, “
Mechanics of Three-Dimensional Printed Lattices for Biomedical Devices
,”
ASME J. Mech. Des.
,
141
(
3
), p.
031703
.
15.
Giannoudis
,
P. V.
,
Einhorn
,
T. A.
, and
Marsh
,
D.
,
2007
, “
Fracture Healing: The Diamond Concept
,”
Injury
,
38
(
Suppl. 4
), pp.
S3
S6
.
16.
Almeida
,
H. A.
, and
Bártolo
,
P. J.
,
2014
, “
Design of Tissue Engineering Scaffolds Based on Hyperbolic Surfaces: Structural Numerical Evaluation
,”
Med. Eng. Phys.
,
36
(
8
), pp.
1033
1040
.
17.
Wieding
,
J.
,
Wolf
,
A.
, and
Bader
,
R.
,
2014
, “
Numerical Optimization of Open-Porous Bone Scaffold Structures to Match the Elastic Properties of Human Cortical Bone
,”
J. Mech. Behav. Biomed. Mater.
,
37
(
21
), pp.
56
68
.
18.
Dias
,
M. R.
,
Guedes
,
J. M.
,
Flanagan
,
C. L.
,
Hollister
,
S. J.
, and
Fernandes
,
P. R.
,
2014
, “
Optimization of Scaffold Design for Bone Tissue Engineering: A Computational and Experimental Study
,”
Med. Eng. Phys.
,
36
(
4
), pp.
448
457
.
19.
Boccaccio
,
A.
,
Fiorentino
,
M.
,
Uva
,
A. E.
,
Laghetti
,
L. N.
, and
Monno
,
G.
,
2018
, “
Rhombicuboctahedron Unit Cell Based Scaffolds for Bone Regeneration: Geometry Optimization With a Mechanobiology Driven Algorithm
,”
Mater. Sci. Eng. C
,
83
(
March 2017
), pp.
51
66
.
20.
Mohammed
,
M. I.
, and
Gibson
,
I.
,
2018
, “
Design of Three-Dimensional, Triply Periodic Unit Cell Scaffold Structures for Additive Manufacturing
,”
ASME J. Mech. Des.
,
140
(
7
), p.
071701
.
21.
Egan
,
P. F.
,
Ferguson
,
S. J.
, and
Shea
,
K.
,
2017
, “
Design of Hierarchical Three-Dimensional Printed Scaffolds Considering Mechanical and Biological Factors for Bone Tissue Engineering
,”
ASME J. Mech. Des.
,
139
(
6
), p.
061401
.
22.
Roberge
,
J.
, and
Norato
,
J.
,
2018
, “
Computational Design of Curvilinear Bone Scaffolds Fabricated Via Direct Ink Writing
,”
CAD Comput. Aid. Des.
,
95
(
12
), pp.
1
13
.
23.
Luo
,
D.
,
Rong
,
Q.
, and
Chen
,
Q.
,
2017
, “
Finite-Element Design and Optimization of a Three-Dimensional Tetrahedral Porous Titanium Scaffold for the Reconstruction of Mandibular Defects
,”
Med. Eng. Phys.
,
47
, pp.
176
183
.
24.
Sutradhar
,
A.
,
Paulino
,
G. H.
,
Miller
,
M. J.
, and
Nguyen
,
T. H.
,
2010
, “
Topological Optimization for Designing Patient-Specific Large Craniofacial Segmental Bone Replacements
,”
Proc. Natl. Acad. Sci. U. S. A.
,
107
(
30
), pp.
13222
13227
.
25.
Makowski
,
P.
, and
Kuś
,
W.
,
2016
, “
Optimization of Bone Scaffold Structures Using Experimental and Numerical Data
,”
Acta Mech.
,
227
(
1
), pp.
139
149
.
26.
Byrne
,
D. P.
,
Lacroix
,
D.
,
Planell
,
J. A.
,
Kelly
,
D. J.
, and
Prendergast
,
P. J.
,
2007
, “
Simulation of Tissue Differentiation in a Scaffold as a Function of Porosity, Young’s Modulus and Dissolution Rate: Application of Mechanobiological Models in Tissue Engineering
,”
Biomaterials
,
28
(
36
), pp.
5544
5554
.
27.
Sanz-Herrera
,
J. A.
,
Garcia-Aznar
,
J. M.
, and
Doblare
,
M.
,
2008
, “
A Mathematical Model for Bone Tissue Regeneration Inside a Specific Type of Scaffold
,”
Biomech. Model. Mechanobiol.
,
7
(
5
), pp.
355
366
.
28.
Bashkuev
,
M.
,
Checa
,
S.
,
Postigo
,
S.
,
Duda
,
G.
, and
Schmidt
,
H.
,
2015
, “
Computational Analyses of Different Intervertebral Cages for Lumbar Spinal Fusion
,”
J. Biomech.
,
48
(
12
), pp.
3274
3282
.
29.
Pobloth
,
A. M.
,
Checa
,
S.
,
Razi
,
H.
,
Petersen
,
A.
,
Weaver
,
J. C.
,
Chmidt-Bleek
,
K.
,
Windolf
,
M.
,
Tatai
,
A. A.
,
Roth
,
C. P.
,
Schaser
,
K. D.
,
Duda
,
G. N.
, and
Schwabe
,
P.
,
2018
, “
Mechanobiologically Optimized 3D Titanium-Mesh Scaffolds Enhance Bone Regeneration in Critical Segmental Defects in Sheep
,”
Sci. Transl. Med.
,
10
(
423
), pp.
1
16
.
30.
Paris
,
M.
,
Götz
,
A.
,
Hettrich
,
I.
,
Bidan
,
C. M.
,
Dunlop
,
J. W.
,
Razi
,
H.
,
Zizak
,
I.
,
Hutmacher
,
D. W.
,
Fratzl
,
P.
,
Duda
,
G. N.
,
Wagermaier
,
W.
, and
Cipitria
,
A.
,
2017
, “
Scaffold Curvature-Mediated Novel Biomineralization Process Originates a Continuous Soft Tissue-to-Bone Interface
,”
Acta Biomater.
,
60
, pp.
64
80
.
31.
Norato
,
J. A.
, and
Wagoner Johnson
,
A. J.
,
2011
, “
A Computational and Cellular Solids Approach to the Stiffness-Based Design of Bone Scaffolds
,”
ASME J. Biomech. Eng.
,
133
(
9
), p.
091003
.
32.
Lan Levengood
,
S. K.
,
Polak
,
S. J.
,
Wheeler
,
M. B.
,
Maki
,
A. J.
,
Clark
,
S. G.
,
Jamison
,
R. D.
, and
Wagoner Johnson
,
A. J.
,
2010
, “
Multiscale Osteointegration as a New Paradigm for the Design of Calcium Phosphate Scaffolds for Bone Regeneration
,”
Biomaterials
,
31
(
13
), pp.
3552
3563
.
33.
Andreassen
,
E.
, and
Andreasen
,
C. S.
,
2014
, “
How to Determine Composite Material Properties Using Numerical Homogenization
,”
Comput. Mater. Sci.
,
83
, pp.
488
495
.
34.
Dong
,
G.
,
Tang
,
Y.
, and
Zhao
,
Y. F.
,
2019
, “
A 149 Line Homogenization Code for Three-Dimensional Cellular Materials Written in MATLAB
,”
ASME J. Eng. Mater. Technol.
,
141
(
1
), p.
011005
.
35.
Terada
,
K.
,
Ito
,
T.
, and
Kikuchi
,
N.
,
1998
, “
Characterization of the Mechanical Behaviors of Solid-Fluid Mixture by the Homogenization Method
,”
Comput. Methods Appl. Mech. Eng.
,
153
(
3–4
), pp.
223
257
.
36.
Norato
,
J.
,
Haber
,
R.
,
Tortorelli
,
D.
, and
Bendsøe
,
M. P.
,
2004
, “
A Geometry Projection Method for Shape Optimization
,”
Int. J. Numer. Methods Eng.
,
60
(
14
), pp.
2289
2312
.
37.
Norato
,
J. A.
,
Bell
,
B. K.
, and
Tortorelli
,
D. A.
,
2015
, “
A Geometry Projection Method for Continuum-Based Topology Optimization With Discrete Elements
,”
Comput. Methods Appl. Mech. Eng.
,
293
(
4
), pp.
306
327
.
38.
Smith
,
H.
, and
Norato
,
J. A.
,
2020
, “
A MATLAB Code for Topology Optimization Using the Geometry Projection Method
,”
Struct. Multidiscipl. Optim.
,
62
, pp.
1579
1594
.
39.
Kazemi
,
H.
,
Vaziri
,
A.
, and
Norato
,
J. A.
,
2018
, “
Topology Optimization of Structures Made of Discrete Geometric Components With Different Materials
,”
ASME J. Mech. Des.
,
140
(
11
), p.
111401
.
40.
Hoelzle
,
D. J.
,
Alleyne
,
A. G.
, and
Wagoner Johnson
,
A. J.
,
2008
, “
Micro-Robotic Deposition Guidelines by a Design of Experiments Approach to Maximize Fabrication Reliability for the Bone Scaffold Application
,”
Acta Biomater.
,
4
(
4
), pp.
897
912
.
41.
Jacobs
,
C. R.
,
1994
, “
Numerical Simulation of Bone Adaptation to Mechanical Loading
,”
PhD thesis
,
Stanford University
,
Stanford, CA
.
42.
Doblaré
,
M.
, and
García
,
J.
,
2001
, “
Application of an Anisotropic Bone-Remodelling Model Based on a Damage-Repair Theory to the Analysis of the Proximal Femur Before and After Total Hip Replacement
,”
J. Biomech.
,
34
(
9
), pp.
1157
1170
.
43.
Fyhrie
,
D. P.
, and
Carter
,
D. R.
,
1986
, “
A Unifying Principle Relating Stress to Trabecular Bone Morphology
,”
J. Orthop. Res.
,
4
(
3
), pp.
304
317
.
44.
Martin
,
R. B.
,
1984
, “
Porosity and Specific Surface of Bone
,”
J. Crit. Rev. Biomed. Eng.
,
10
(
3
), pp.
179
222
.
45.
Greenwald
,
A. S.
,
Boden
,
S. D.
,
Goldberg
,
V. M.
,
Khan
,
Y.
,
Laurencin
,
C. T.
, and
Rosier
,
R. N.
,
2001
, “
Bone-Graft Substitutes: Facts, Fictions, and Applications
,”
J. Bone Joint Surg.
46.
Beaupré
,
G. S.
,
Orr
,
T. E.
, and
Carter
,
D. R.
,
1990
, “
An Approach for Time-Dependent Bone Modeling and Remodeling-Application: A Preliminary Remodeling Simulation
,”
J. Orthop. Res.
,
8
(
5
), pp.
663
670
.
47.
Beaupré
,
G. S.
,
Orr
,
T. E.
, and
Carter
,
D. R.
,
1990
, “
An Approach for Time-Dependent Bone Modeling and Remodeling- Theoretical Development
,”
J. Orthop. Res.
,
8
(
5
), pp.
651
661
.
48.
Liu
,
Y.
,
Zheng
,
G.
,
Letov
,
N.
, and
Zhao
,
Y. F.
,
2021
, “
A Survey of Modeling and Optimization Methods for Multi-Scale Heterogeneous Lattice Structures
,”
ASME J. Mech. Des.
,
143
(
4
), p.
040803
.
49.
Zhang
,
S.
,
Norato
,
J. A.
,
Gain
,
A. L.
, and
Lyu
,
N.
,
2016
, “
A Geometry Projection Method for the Topology Optimization of Plate Structures
,”
Struct. Multidiscipl. Optim.
,
54
(
5
), pp.
1173
1190
.
50.
Hollister
,
S. J.
, and
Kikuchi
,
N.
,
1992
, “
A Comparison of Homogenization and Standard Mechanics Analyses for Periodic Porous Composites
,”
Comput. Mech.
,
10
(
2
), pp.
73
95
.
51.
Zhu
,
T.
,
Cui
,
Y.
,
Zhang
,
M.
,
Zhao
,
D.
,
Liu
,
G.
, and
Ding
,
J.
,
2020
, “
Engineered Three-Dimensional Scaffolds for Enhanced Bone Regeneration in Osteonecrosis
,”
Bioact. Mater.
,
5
(
3
), pp.
584
601
.
52.
Weiser
,
A.
, and
Zarantonello
,
S. E.
,
1988
, “
A Note on Piecewise Linear and Multilinear Table Interpolation in Many Dimensions
,”
Math. Comput.
,
50
, pp.
181
189
.
53.
Viana
,
F. A.
,
Gogu
,
C.
, and
Haftka
,
R. T.
,
2010
, “
Making the Most Out of Surrogate Models: Tricks of the Trade
,”
Proc. ASME Des. Eng. Tech. Conf.
,
1
(PARTS A AND B), pp.
587
598
.
54.
Afzal
,
A.
,
Kim
,
K. Y.
, and
Seo
,
J. W.
,
2017
, “
Effects of Latin Hypercube Sampling on Surrogate Modeling and Optimization
,”
Int. J. Fluid Mach. Syst.
,
10
(
3
), pp.
240
253
.
55.
Amini
,
A. R.
,
Laurencin
,
C. T.
, and
Nukavarapu
,
S. P.
,
2012
, “
Bone Tissue Engineering: Recent Advances and Challenges
,”
Crit. Rev. Biomed. Eng.
,
40
(
5
), pp.
363
408
.
56.
Karageorgiou
,
V.
, and
Kaplan
,
D.
,
2005
, “
Porosity of 3d Biomaterial Scaffolds and Osteogenesis
,”
Biomaterials
,
26
(
27
), pp.
5474
5491
.
57.
Razi
,
T.
,
Niknami
,
M.
, and
Alavi Ghazani
,
F.
,
2014
, “
Relationship Between Hounsfield Unit in CT Scan and Gray Scale in CBCT.
,”
J. Dental Res. Dental Clin. Dental Prospects
,
8
(
2
), pp.
107
110
.
58.
Kaneko
,
T. S.
,
Pejcic
,
M. R.
,
Tehranzadeh
,
J.
, and
Keyak
,
J. H.
,
2003
, “
Relationships Between Material Properties and CT Scan Data of Cortical Bone With and Without Metastatic Lesions
,”
Med. Eng. Phys.
,
25
(
6
), pp.
445
454
.
59.
Schileo
,
E.
,
Dall’Ara
,
E.
,
Taddei
,
F.
,
Malandrino
,
A.
,
Schotkamp
,
T.
,
Baleani
,
M.
, and
Viceconti
,
M.
,
2008
, “
An Accurate Estimation of Bone Density Improves the Accuracy of Subject-Specific Finite Element Models
,”
J. Biomech.
,
41
(
11
), pp.
2483
2491
.
60.
Ohs
,
N.
,
Collins
,
C. J.
, and
Atkins
,
P. R.
,
2020
, “
Validation of HR-pQCT Against Micro-CT for Morphometric and Biomechanical Analyses: A Review
,”
Bone Reports
,
13
(
Suppl. 5
), p.
100711
.
61.
Metcalf
,
L. M.
,
Dall’Ara
,
E.
,
Paggiosi
,
M. A.
,
Rochester
,
J. R.
,
Vilayphiou
,
N.
,
Kemp
,
G. J.
, and
McCloskey
,
E. V.
,
2018
, “
Validation of Calcaneus Trabecular Microstructure Measurements by HR-pQCT
,”
Bone
,
106
(
6
), pp.
69
77
.
62.
Burghardt
,
A. J.
,
Kazakia
,
G. J.
, and
Majumdar
,
S.
,
2007
, “
A Local Adaptive Threshold Strategy for High Resolution Peripheral Quantitative Computed Tomography of Trabecular Bone
,”
Ann. Biomed. Eng.
,
35
(
10
), pp.
1678
1686
.
63.
Nazarian
,
A.
,
Snyder
,
B. D.
,
Zurakowski
,
D.
, and
Müller
,
R.
,
2008
, “
Quantitative Micro-Computed Tomography: A Non-Invasive Method to Assess Equivalent Bone Mineral Density
,”
Bone
,
43
(
2
), pp.
302
311
.
64.
Dassault Systems
,
2019
, “
Abaqus 3DEXPERIENCE R2019x
.”
65.
Gushue
,
D. L.
,
Houck
,
J.
, and
Lerner
,
A. L.
,
2005
, “
Rabbit Knee Joint Biomechanics: Motion Analysis and Modeling of Forces During Hopping
,”
J. Orthop. Res.
,
23
(
4
), pp.
735
742
.
66.
Bouhlel
,
M. A.
,
Hwang
,
J. T.
,
Bartoli
,
N.
,
Lafage
,
R.
,
Morlier
,
J.
, and
Martins
,
J. R. R. A.
,
2019
, “
A Python Surrogate Modeling Framework With Derivatives
,”
Adv. Eng. Softw.
,
135
(
5
), p.
102662
.
67.
Virtanen
,
P.
,
Gommers
,
R.
,
Oliphant
,
T. E.
,
Haberland
,
M.
,
Reddy
,
T.
,
Cournapeau
,
D.
,
Burovski
,
E.
,
Peterson
,
P.
,
Weckesser
,
W.
,
Bright
,
J.
,
van der Walt
,
S. J.
,
Brett
,
M.
,
Wilson
,
J.
,
Millman
,
K. J.
,
Mayorov
,
N.
,
Nelson
,
A. R. J.
,
Jones
,
E.
,
Kern
,
R.
,
Larson
,
E.
,
Carey
,
C. J.
,
Polat
,
İ.
,
Feng
,
Y.
,
Moore
,
E. W.
,
VanderPlas
,
J.
,
Laxalde
,
D.
,
Perktold
,
J.
,
Cimrman
,
R.
,
Henriksen
,
I.
,
Quintero
,
E. A.
,
Harris
,
C. R.
,
Archibald
,
A. M.
,
Ribeiro
,
A. H.
,
Pedregosa
,
F.
,
van Mulbregt
,
P.
, and
SciPy 1.0 Contributors
,
2020
, “
SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python
,”
Nat. Methods
,
17
(
3
), pp.
261
272
.
68.
Polak
,
S. J.
,
Levengood
,
S. K.
,
Wheeler
,
M. B.
,
Maki
,
A. J.
,
Clark
,
S. G.
, and
Johnson
,
A. J.
,
2011
, “
Analysis of the Roles of Microporosity and BMP-2 on Multiple Measures of Bone Regeneration and Healing in Calcium Phosphate Scaffolds
,”
Acta Biomater.
,
7
(
4
), pp.
1760
1771
.
69.
Lan Levengood
,
S. K.
,
Polak
,
S. J.
,
Poellmann
,
M. J.
,
Hoelzle
,
D. J.
,
Maki
,
A. J.
,
Clark
,
S. G.
,
Wheeler
,
M. B.
, and
Johnson
,
A. J.
,
2010
, “
The Effect of BMP-2 on Micro- and Macroscale Osteointegration of Biphasic Calcium Phosphate Scaffolds With Multiscale Porosity
,”
Acta Biomater.
,
6
(
8
), pp.
3283
3291
.