Medical Engineering & Physics 36 (2014) 448–457

Contents lists available at ScienceDirect

Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy

Optimization of scaffold design for bone tissue engineering: A computational and experimental study Marta R. Dias a , José M. Guedes a , Colleen L. Flanagan b , Scott J. Hollister b , Paulo R. Fernandes a,∗ a b

IDMEC, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais, 1, 1049-001 Lisboa, Portugal Scaffold Tissue Engineering Group, University of Michigan, 3412 GGB, 2350 Hayward, Ann Arbor, MI 48109-2125, USA

a r t i c l e

i n f o

Article history: Received 21 May 2013 Received in revised form 2 December 2013 Accepted 6 February 2014 Keywords: Topology optimization Homogenization Bone tissue engineering Selective laser sintering

a b s t r a c t In bone tissue engineering, the scaffold has not only to allow the diffusion of cells, nutrients and oxygen but also provide adequate mechanical support. One way to ensure the scaffold has the right properties is to use computational tools to design such a scaffold coupled with additive manufacturing to build the scaffolds to the resulting optimized design specifications. In this study a topology optimization algorithm is proposed as a technique to design scaffolds that meet specific requirements for mass transport and mechanical load bearing. Several micro-structures obtained computationally are presented. Designed scaffolds were then built using selective laser sintering and the actual features of the fabricated scaffolds were measured and compared to the designed values. It was possible to obtain scaffolds with an internal geometry that reasonably matched the computational design (within 14% of porosity target, 40% for strut size and 55% for throat size in the building direction and 15% for strut size and 17% for throat size perpendicular to the building direction). These results support the use of these kind of computational algorithms to design optimized scaffolds with specific target properties and confirm the value of these techniques for bone tissue engineering. © 2014 IPEM. Published by Elsevier Ltd. All rights reserved.

1. Introduction Scaffolds for bone tissue engineering must be highly permeable in order to promote cell proliferation and differentiation as well as allowing oxygen, nutrient and metabolic waste diffusion, but they have also to provide the necessary mechanical support [1–3], which will depend upon anatomic location. Assuming bone is adapted to its load bearing function, the elastic constants of normal healthy bone can provide a design target for a specific anatomic application. Under such assumption, there will be circumstances where the elastic properties should be maximized, but other applications where only a minimum of stiffness may be required, as in completely contained defects [4]. In addition, permeability has been demonstrated in many cases to be associated with higher bone regeneration [5–8]. Indeed, permeability is important for bone regeneration not only because higher values improve bone ingrowth but also because inadequate values may induce the

∗ Corresponding author at: IDMEC-IST, Av. Rovisco Pais, 1, 1049-001 Lisboa, Portugal. Tel.: +351 218417925; fax: +351 218417915. E-mail addresses: [email protected] (M.R. Dias), [email protected] (J.M. Guedes), fl[email protected] (C.L. Flanagan), [email protected] (S.J. Hollister), [email protected] (P.R. Fernandes). http://dx.doi.org/10.1016/j.medengphy.2014.02.010 1350-4533/© 2014 IPEM. Published by Elsevier Ltd. All rights reserved.

formation of cartilaginous tissue instead of bone [9,10]. Several studies on bone scaffold have focused on permeability as a key parameter [11,12], including a number of authors who have been focused on the development of computational tools, as well as measuring systems for permeability analysis [11,13–16]. The objective of this work is to develop an optimization tool able to design scaffolds that can achieve a target performance with respect to stiffness and permeability, having in mind the goal of designing scaffolds for specific clinical applications. The problem consists of a material distribution problem based on topology optimization of structures [17]. The idea of designing the scaffold microstructure based on topology optimization has also been explored by other researchers. Some of the first studies were presented by Hollister and co-workers [6,18]. It was demonstrated that tissue regeneration depends on scaffold architecture parameters like porosity and permeability and consequently the control of scaffold architecture may help explain the relationship of scaffold architecture to biological behavior. Although it was not for a tissue engineering application, Guest and Prévost have also proposed a work on topology optimization to design periodic porous structures with maximum effective permeability and prescribed flow properties [19]. More recently, Li’s research group has presented studies on scaffold design based on topology optimization [20,21].

M.R. Dias et al. / Medical Engineering & Physics 36 (2014) 448–457

In the present study a topology optimization approach is presented as a tool to design bone scaffolds for given mechanical conditions (strain field) and/or required permeability. Several computational results, concerning different optimization scenarios, will be presented. Some of the resulting designs were built from an implantable biomaterial, polycaprolactone-4%hydroxyapatite (PCL-4%HA), using a selective laser sintering (SLS) system [22] and the actual features of the scaffolds, including the pore interconnection sizes (here called throats for simplification), strut sizes and porosity, were assessed experimentally in comparison to the optimally designed features. The paper represents a departure from previous works since here a complete study including scaffold optimization for elasticity and permeability, SLS fabrication and experimental validation is done. The study of this chain from design to manufacture demonstrates the unique challenges inherent in optimizing scaffold architecture and subsequently fabricating the scaffold to replicate the original complex, optimized design. 2. Materials and methods 2.1. Scaffolds design – topology optimization Topology optimization consists of defining pore and material regions within a given design domain [17]. Here, the design domain is a unit cell which repeated periodically in the three-dimensional space defines the scaffold microstructure, assumed as a periodic porous media (Fig. 1). The quantity of material allowed to distribute within the unit cell defines the volume fraction (or inversely the porosity) of the scaffold. In this work it is assumed scaffolds are homogeneous with periodic structure. However, the scaffold could be heterogeneous and designed “region by region” [23]. The permeability and elastic properties of this periodic media are computed using an asymptotic homogenization method as described for instance by Guedes and Kikuchi [24]. To compute the permeability properties, the problem of a fluid flowing through a porous media is described by a 2nd order differential equation which is obtained by a combination of the Darcy Law and the balance of steady state flow [14]. Its integral form is described by Eq. (1), where K is the second order tensor of the medium permeability coefficients (m2 ), P is the hydraulic pressure (Pa),  is the fluid viscosity (Pa s), Ф is an admissible arbitrary smooth weight function (see for example [25]), f is the quantity of fluid being removed or generated by volume (s−1 ) and q is the Darcy flux on the boundary  q (m/s):



Kij ˝

∂P d d˝ − ∂xj ∂xi





fd˝ − ˝

∀ admissible

qd = 0, q

(1)

Other authors have derived Darcy’s Law from homogenization [19,26,27], with the difference that these have utilized Stokes flow

as the basis for expanding the homogenization description whereas here Darcy’s law is directly used to expand the description. The formulation in Eq. (1), however, provides for easier computation of sensitivity derivatives and is thus easier to implement in the topology optimization algorithm. For elastic properties let us consider the elasticity problem in its integral form (see for example [25]), where E is the fourth order tensor of elastic coefficients (Pa), u is the displacement field (m), v is the virtual displacement field, b represents the body forces (N/m3 ) and t the traction on the boundary t (N/m2 ).



Eijkm ˝

∂uk ∂vi d˝ − ∂xm ∂xj





bi vi d˝ − ˝

ti vi d,

(2)

i

∀vi admissible

The homogenization method “upscales” the problems described by Eqs. (1) and (2), converting the domain into a homogenized equivalent one, avoiding the complex task of analyzing the details of a periodic material. Using the method described by Guedes and Kikuchi [24], the equivalent permeability and elastic coefficients are obtained, respectively, by: H kim

1 = Y

H Eijkm

=





Kij

∂xm ıjm − ∂yj





y

1

  Y 

Epqrs ¥



dY

∂x¯ km ∂rk ∂sm r ∂ys

(3)

 ∂pi ∂qj

ij

∂x¯ p

∂yq

 d¥

(4)

The  functions in Eq. (3) represent the microstructure pressure perturbations for a unit average pressure gradient in each direction and are the solution of a series of problems in the microstructure:





Kij Y

∂xm ∂ dY = − ∂yj ∂yi

kim Y

∂ dY ∂yi

∀ admissible

(5)

Equivalently, the x¯ functions, in Eq. (4), represent the deformation modes for a unit cell subject to six unit average strains (3 normal and 3 shear), given by the solutions of the following set of problems in the microstructure:



Eijrs ¥

∂ ¯ km ∂vi r = ∂ys ∂yj



Eijkm ¥

∂ vi dY ∂yj

∀ vi admissible

(6)

These problems are solved on the domain of the unit cell by the finite elements method (FEM), using a MATLAB code. The unit cell is discretized in a 20 × 20 × 20 mesh of 8-node brick elements. To define the topology of the unit cell, a continuous “density field” ␳ is associated to each finite element and its value is 1 if the finite element corresponds to a material region and 0 if it is void. This means that each element has a characterizing density value and, based on that value, each element’s mechanical and permeability properties are defined by a power law, where K0 = 1, E0 = 10 and n is chosen as between 4 and 6: k = (1 − )n K 0 E = n E 0

Fig. 1. Homogenization theory. Left, domain of the scaffold ε ; center, detail of the domain; right, unit cell Y; d, characteristic length scale of the microstructure size, represented by y; D, characteristic length scale of the scaffold, represented by x; ε, ratio between d and D.

449

(7)

The final geometry will preferably only have 0 and 1 densities corresponding to material and void elements. In fact, this way of defining the material law corresponds to the SIMP (solid isotropic material with penalization) model [17], which leads the design variable to its extreme values, 0 or 1, once a suitable exponent n is chosen. Material elements will have null permeability and maximum elasticity and void elements the other way around. The homogenized properties are normalized to K0 , the “void permeability” and E0 , the base material elastic properties. While to obtain the effective elasticity, one only needs to multiply the resulting E by the value of the elasticity of the base material (Ebase ), for permeability it is more complex since the permeability for the

450

M.R. Dias et al. / Medical Engineering & Physics 36 (2014) 448–457

void region is not easy to define [14]. Thus, after obtaining the scaffold’s internal geometry by optimization, the final permeability properties values are computed by homogenizing the Stokes equation. Here the homogenized Stokes equation was considered a good approach since it is a onetime calculation after the optimization process and it allows permeability computation based on the real dimensions of the unit cell. The equation can be obtained by solving the microstructure Stokes flow due to unit average pressure gradients (see e.g. [3]):

 Yf

j

∂ki ∂ϕi dYf = ∂yp ∂yp

 ϕj dYf

(8)

1

  Y 



j

Yf

ki dYf

(9)

The permeability function will be defined as the average of the homogenized coefficients on the main directions and the elasticity function will be defined as the elastic strain energy, which depends on homogenized elastic coefficients and on the applied strain field e: Fperm =

H + KH + KH K11 22 33

3

Felast

(10)

1 = eE H e 2

Here KH are the homogenized permeability coefficients given by Eq. (3) and EH are the homogenized elastic coefficients given by Eq. (4). Assuming the material model described above, the optimization problem can be stated based on the aims of the scaffolds application. In this study two problems are defined. In both, the objective function is the inverse of Fperm (given by Eq. (10)), but the applied constraints differ. In the first case, a lower bound is imposed on the elastic strain ∗ energy (Felast ). The constraint value (Felast ) will be defined as a percentage of the maximum elastic strain energy (situation equivalent to a full material unit cell). The topology optimization problem is thus described by Eq. (11), where is the density field: Min

1 Fperm

∗ subject to Felast ≥Felast

(11)

≥0, ≤ 1 Results for two different strain fields, e = (1,1,1,0,0,0) and ∗ e = (0,0,0,1,1,1), with Felast = 0.30, as well as for four different con∗ straint values, Felast = 0.01, 0.20, 0.40, 0.60, with e = (1,1,1,0,0,0), will be presented. In the second case, a lower bound is imposed on volume fraction (VF), in addition to constraints on the elastic strain energy, as before. The topology optimization problem is then defined by: Min

2.2. Data preparation

Yf

And permeability is thus given by: KijH =

consequently the same mechanical properties but different degrees of permeability due to the different throat sizes resulted from the homothetic expansion. The optimization method chosen to solve the proposed topology optimization problem was the method of moving asymptotes (MMA) [28]. In the MMA, in each iteration, a separable and convex sub-problem is generated (based on the gradients of the objective function and constraints of the original problem) and thus dual methods can be used, reducing the computational effort.

1 Fperm

∗ subject to Felast ≥Felast

(12)

The resulting design for the second case was used to manufacture scaffolds with the purpose of studying their actual features in comparison to the prescribed ones. After obtaining the optimized microstructure, we “binarized” the few elements that remain with density between zero and one by assigning 1 or 0 using a threshold value, typically 0.5. The FE mesh with 20 × 20 × 20 elements was converted into a stack of 2D .tiff images, by slicing the 3D mesh using a MATLAB algorithm (Fig. 2). The number of slices and the number of unitary pixels in each of the 2D image were defined in order to homothetically expand the unit cell to three different sizes, 2 mm (2p0), 3 mm (3p0) and 4 mm (4p0). Thus, for the 2p0 unit cell there were 80 .tiff images with 80 × 80 voxels, for the 3p0 unit cell there were 60 .tiff images with 60 × 60 unitary voxels and for the 4p0 unit cell there 80 .tiff images with 80 × 80 voxels. These sets of images were in turn converted to 3D voxel based datasets using an IDL program (Interactive Data LanguageTM , IDL; Research Systems, Inc., Boulder, CO, USA), where the voxel size was defined as 0.025 mm for the 2p0 and 0.05 mm for the 3p0 and 4p0 designs, thus achieving the sizes of 2, 3 and 4 mm for the unit cell (Fig. 2). The 3D cylindrical scaffolds were then designed based on assembling repetitions of the unit cell, using another customized IDL program (Fig. 2). The 3-dimensional image based dataset was converted to a surface based geometry STL file using an IDL program (Fig. 2), based on the Marching Cubes algorithm, where the voxels are split into tetrahedrals based on an isodensity surface and each triangle vertex location normal vector projected in the STL file is “relaxed” by averaging it with surrounding vertices. This is done iteratively using a weighting function, where the weighting function and number of iterations was fixed. 2.3. Scaffold manufacturing The scaffolds were built from poly-␧-caprolactone4%hydroxyapatite (Polysciences, Warrington, PA; Plasma Biotal Limited, UK), using a selective laser sintering system (EOS Formiga P 100), based on the STL files. The .STL files were prepared and sliced for rapid prototyping using RPTools software based on a layer thickness of 0.08 mm. The polycaprolactone was cryogenically milled, such that it had a mean particle size range of 40–50 ␮m with a maximum particle size of 125 ␮m, and then blended with hydroxyapatite. The scaffolds were manufactured in a nitrogen environment, with a bed temperature range of 52–55 ◦ C and laser power and speed settings of 4 W and 1500 mm/s, respectively.

VF≥VF ∗ ≥0, ≤ 1 The result is an optimized material distribution, which defines the geometry of the unit cell. After, three different unit cell sizes were chosen and the resulting geometry was homothetically expanded accordingly. This way, it was possible to obtain three scaffolds, all with the same geometry, the same volume fraction and

2.4. Scaffold feature measurement To measure the actual features of the manufactured scaffolds, Hounsfield calibrated micro-CT scans were generated using a high-resolution micro-CT scanner (eXplore Locus, GE Healthcare, Milwaukee, WI, kVp = 80, uA = 450, integration time = 2000 ms) with a voxel resolution of 20.7 ␮m.

M.R. Dias et al. / Medical Engineering & Physics 36 (2014) 448–457

451

Fig. 2. From design to manufacture: initial and optimized FE meshes (dark gray – void, light gray – material), set of .tiff images, STL image of the scaffold with details of 4 unit cells (light gray – material), and manufactured scaffolds (white – material).

For the porosity measurement a region of interest (ROI) that encompassed a unit cell volume was selected (an average of 3 unit cells per scaffold) and, using a fixed threshold level of −600 HU to identify the material, the solid volume fraction and consequently the porosity values were calculated using MicroView software (GE Medical Systems, Toronto, Canada). An average of three scaffolds was measured to get an average of 8 measurements for each design. For the throat measurements, one cross section perpendicular to each direction was chosen from the interior of the scaffold to capture one of the smallest interconnections of the scaffold (Fig. 3). The equivalent was done for the strut measurements, choosing a cross section that captured one of the smallest strut diameters from the interior of the scaffold (Fig. 3). From each cross section, a minimum of two struts and throats were measured in both the horizontal and vertical directions. An average of three scaffolds was measured for each design. At the end, each feature in each direction (e.g. strut horizontal size in z-direction) will be given by the average of five measurements. Moreover, because the cylindrical scaffolds were

built with their long axis aligned with the z-direction of the SLS machine, the features will be described in that direction and on the perpendicular plane to that direction, x/y plane. Thus the measurements taken at the planes perpendicular to y and x-direction were averaged, as it was not possible to uniquely discern the x and ydirections as they were identical. All the measurements were done using MicroView software. Also, a ROI corresponding to just one layer, encompassing a unit cell, either centered on one of the smallest struts or smallest throats, was selected and its volume fraction was used to compute the struts and throats areas respectively. 3. Results 3.1. Optimal scaffold geometries First, the objective function was defined as Fperm (Eq. (10)) and elastic strain energy was constrained to be greater than a lower bound to ensure material continuity and thus the ability to

Fig. 3. Schematic of the measurement system: in z-direction (left) with details of one of the smallest throats, in y-direction (center), in x-direction (right) with details of one of the smallest struts, white – material, black – void.

452

M.R. Dias et al. / Medical Engineering & Physics 36 (2014) 448–457

Fig. 4. Resulting unit cells for the first case, for two strain fields: Left, e = (1,1,1,0,0,0), an entire unit cell and a sliced unit cell; right, e = (0,0,0,1,1,1), an entire unit cell and a ∗ representation of just the material, with Felast = 0.30 (white – material, black – void).

actually manufacture the scaffolds. Two conditions were investigated, the application of different strain fields and different degrees of minimum elastic strain energy. The resulting structures present different geometries, depending on the strain field (Fig. 4). If a normal field is applied, the material will be distributed along the three main directions, originating solid struts on those directions. On the other hand, if a shear field is applied, the material will aggregate mainly on a diagonal direction. An interesting feature of these two types of structures is that regardless of the direction of the solid struts, the intrinsic permeability is guaranteed as a consequence of the maximization of Fperm . To study the influence of different values for the elastic strain energy, a normal field was applied and four different values were ∗ imposed as constraints (Felast = 0.01, 0.20, 0.40, 0.60). The resulting structures presented similar geometries, with solid struts in the three main directions, but as the constraint value on elastic strain energy increased, the sizes of those struts and the solid volume fraction also increased, from 7.54% to 88.12% (Fig. 5). Again, these structures have permeability in all directions (from 0.86*K0 to 0.047*K0 ) and this is probably related to the way the optimization problem is defined, as permeability is the objective function. For the structures of the first scenario the porosity is the result of the imposed stiffness. Because the study of structures with equal porosity but diverse permeability is important to assess the influence of permeability on bone growth, the next scenario studied structures with maximized permeability, but with constrained volume fraction as well as elastic strain energy. For a volume fraction

constraint of VF = 45.9%, a strain field of e = (1, 1, 1, 0, 0, 0) and ∗ a constraint of Felast = 0.30, the resulting micro-structure is presented in Fig. 2b. Then, the unit cell was homothetically expanded to a size of 2 mm 2p0 design), 3 mm (3p0 design) and 4 mm (4p0 design). The minimum unit cell size of 2 mm was necessary in order to resolve features when manufacturing on the SLS machine. The final result is thus three different scaffolds with the same geometry (Fig. 2left), the same volume fraction (VF = 45.9%) and consequently the same mechanical properties (66.45 MPa), but with different prescribed permeability values (7.3828 × 10−09 , 1.6611 × 10−08 , 2.9530 × 10−08 m2 ). 3.2. Scaffold actual features Fig. 2 presents the main steps from design to manufacture: the voxel FE mesh of the optimized unit cell, the stack of .tiff images, the STL file generated by isodensity surfaces using Marching Cubes, and the manufactured scaffolds. The actual features of the manufactured scaffolds (porosity, throat and strut sizes) were measured based on the micro-CT images and compared with the designed features, as well as with the STL images. The measured features are summarized in Table A1, on the Appendix, and compared in Figs. 6–11. Porosity is slightly higher on the manufactured scaffolds, compared to the STL and designed values (Fig. 6). The actual porosity of the 2p0 design is 13.87% higher than the prescribed value, while for 3p0 and 4p0, it is within 8% of the designed values (Table 1).

Fig. 5. Optimized scaffolds resulting from constraining the elastic energy to 4 different values: from left to right – 1%, 20%, 40% and 60% of the maximum elastic energy, with resulting permeability of K = 0.86 × K0 , 0.30 × K0 , 0.13 × K0 , 0.047 × K0 (white – material, black – void).

M.R. Dias et al. / Medical Engineering & Physics 36 (2014) 448–457

453

This may be due to presence of small defects and the material microporosity, which was not accounted for on the computational calculations (see Fig. 9). The struts and throat sizes of the manufactured scaffolds are slightly different from the designed ones. The struts sizes in z-direction are smaller for the manufactured scaffolds, both horizontally and vertically (Fig. 7), and the throats sizes in z-direction, are slightly bigger for the manufactured scaffolds (both horizontally and vertically), when compared to the designed ones (Fig. 8). In the x/y-direction, the struts are in general slightly smaller when measured horizontally and slightly bigger then measured vertically (Fig. 7). Conversely, the throats are bigger horizontally and smaller when measured vertically (Fig. 8). This means there is Fig. 6. Porosity of scaffolds – designed, STL and manufactured values.

Fig. 7. Minimum struts sizes in z and x/y-direction – designed, STL and actual measurements.

Fig. 8. Minimum throats sizes in z and x/y-direction – designed, STL and actual measurements.

Fig. 9. Thresholded ␮CT images of the struts (top) and throats (bottom), in z and x/y-direction of 2p0, 3p0 and 4p0 (white-material, black–void).

454

M.R. Dias et al. / Medical Engineering & Physics 36 (2014) 448–457

Fig. 10. Strut areas in z and x/y-direction of the manufactured scaffolds, based on microCT images.

Fig. 11. Throat areas in z and x/y-direction of the manufactured scaffolds, based on microCT images.

Table 1 Differences (in %) between designed and actual features (porosity, struts and throats areas). Diff. in porosity (%)

2p0 3p0 4p0

13.87 ± 3.05 4.76 ± 2.13 7.09 ± 3.18

Diff. in strut area (%)

Diff. in throat area (%)

In z-direction

In x/y-direction

In z-direction

In x/y-direction

39.40 ± 4.63 29.02 ± 5.70 21.45 ± 4.88

14.404 ± 5.13 15.15 ± 7.90 0.26 ± 2.92

54.98 ± 13.54 4.02 ± 4.54 19.91 ± 7.56

16.18 ± 2.37 10.19 ± 2.32 10.68 ± 3.46

some distortion on struts and throats in the x/y-direction (compare Fig. 9 with Fig. 2). As shown in Fig. 9, the struts maintain a circular shape in the z-direction, as the designed structures, while these look more vertically distorted in the x/y-direction. The throats have a more diamond-type shape, losing some of the definitions of the crosslike shape it had on design structure (see Fig. 2), which is more evident in the z-direction. Also, it might be important to notice that this discrepancy in shape is more evident for the smallest unit cell, probably due to the fact that this had features closer to the limit of the machine’s resolution. The struts areas are smaller in z-direction (within 40% of the designed value), but in x/y-direction these are slightly bigger (less than 15.2%) than the designed values (Fig. 10 and Table 1). The actual throat areas are slightly bigger in any direction, within 20% of the designed values, except for the 2p0 design (Fig. 11 and Table 1). 4. Discussion Results for the first scenario show that for different applied strain, different optimized structures are obtained (see Fig. 4). In both, the material was distributed in a way that leads to the best response to the field applied. In fact, the possibility of obtaining

scaffolds for different mechanical conditions is one of the advantages of the methodology proposed in this work, which may be important when designing scaffolds for different clinical applications. For instance, the mechanical environment in a non-loaded area (to fill a bone cyst or to repair a cranial defect) is rather different from what is observed on a segmental diaphyseal defect reconstruction [4]. Thus, bone substitutes used in those situations may have different micro-structures in order to best respond to the mechanical requirements. Another aspect is that, because the permeability properties were being maximized, it could be assured that these scaffolds have the best diffusion properties, regardless of the strain field applied. As observed in Fig. 5, as one increases the minimum value of the prescribed elastic energy, the amount of material also increases, since there is not a constraint on volume fraction. Accordingly, the permeability decreases as the elasticity properties increase. This is an example of how this approach can achieve different solutions based on the mechanical requirements. In the second scenario the problem was formulated to achieve scaffolds with different permeability values but the same volume fraction and mechanical properties, in order to use it in future studies on the influence of permeability on bone growth, similar to the one presented by Mitsak [7], although in that case

M.R. Dias et al. / Medical Engineering & Physics 36 (2014) 448–457

changes in permeability were achieved by changing volume fraction but keeping unit cell topology the same. The resulting structure was similar to the one obtained for the first scenario, but with a controlled prescribed volume fraction. With the homothetic expansion of the unit cell, three different sizes were achieved (2 mm, 3 mm and 4 mm) and consequently each design will have different throat sizes. Therefore, it was possible to obtain three scaffolds, with the same internal geometry, the same volume fraction (45.9%) and consequently the same mechanical properties (66.45 MPa), but with different values of permeability (7.4 × 10−9 , 1.7 × 10−8 , 3.0 × 10−8 m2 ). It was thus verified that this approach allows the design of scaffolds for a required value of permeability, just by expanding the unit cell and thus the throat sizes. An interesting aspect of the present work is the constraint on stiffness being introduced through the elastic strain energy. Different approaches can be considered such as direct constraints on Young’s modulus along one direction or the value of elastic constants of the elasticity tensor, alone or combined to satisfy some material symmetry requirements. However, the elastic strain energy approach is more generic because it provides the stiffest structure for the considered strain field. This way it is possible to consider the actual strain field existing on bone and design the scaffold in order to have the best (in energy sense) material properties. This is a natural approach for multiscale analysis as proposed in Coelho et al. [29]. The next question is how the manufactured scaffolds match the designed structures. The complexity of the optimized scaffold pore architecture requires the use of additive manufacturing/3D printing techniques. Each of these techniques have inherent feature limitations due to the form in which the material is used and the physics of how the material is processed. In laser sintering, the method used in this work, there is a minimum feature size in the z-direction (or build direction) based on the layer thickness and particle size distribution that may be milled and spread in the machine. In addition, the laser beam used for sintering has a finite beam diameter (and an effective beam diameter based on the laser power, speed and interaction with the build material), and therefore the feature size in the x–y-directions cannot be smaller than this effective beam width, which for our system and laser parameter settings is 700 ␮m. Some of the features of the scaffold, such as porosity, strut sizes, and throat sizes were measured based on the micro-CT images of the scaffolds.

455

The differences between designed and STL features (Figs. 6–8) were very small and consistent. These are most likely due to the conversion of the FE mesh to a surface based file (see Fig. 2): there are no errors associated from the conversion of the FE mesh to .tiff images since this is done in a directly proportional way, but the smoothing during the conversion to an STL file has a small error associated. However, these errors are not too significant which is confirmed by the small differences on porosity values, struts sizes and throat sizes and as it is possible to observe in Figs. 6–8. However, when comparing the experimental data derived from the micro-CT analysis with the designed values, it was possible to observe some differences not only on the values of porosity and sizes of the struts and throats, but also on the shapes of these two features. The strut areas in the z-direction were smaller than the designed ones (up to 40% of the prescribed values), but the struts areas in the other directions were slightly bigger, up to 15.2%, than the designed values (Fig. 10). These discrepancies are probably related to the differences in shape along the different directions observed in Fig. 9: while in z-direction the struts maintained the circular shape, which associated with smaller struts sizes leads to smaller areas, in the x/y-direction there is a distortion of the struts and thus, although the struts were generally smaller when measured horizontally, the overall areas were slightly bigger. In terms of actual throat areas, these were slightly bigger in all directions, within 20% of the designed values, with exception of the 2p0 design (see Fig. 11). An interesting observation is that in general the 2p0 design presented the highest discrepancies, which are probably due to their features being closer to the SLS resolution that may directly influence the manufacturing accuracy and due to the fact that their small interconnections may have been compromised in the depowdering process. The highest values observed for the porosity of the manufactured scaffolds are probably due to the microporosity of the material, also called “manufacturing-induced porosity”, which contributes to the overall porosity of the scaffold. This effect has been observed by other authors [30,31] and it has been associated to the incomplete coalescence of polymer particles during sintering, melting and resolidification. However, microporosity inherent to SLS is quite low. In a previous study the microporosity for solid PCL cylindrical specimens was less than 1% and for macroporous 3D specimens it was less than 4% [31]. In fact, the differences between designed and actual features in scaffolds built by SLS have been reported in other studies, e.g.

Fig. 12. FE meshes (top), struts (center) and throats (bottom), in z and x/y-direction of 2p0 (left), 3p0 (center) and 4p0 (right), white – material, black – void.

456

M.R. Dias et al. / Medical Engineering & Physics 36 (2014) 448–457

Fig. 13. Permeability (left) and mechanical properties (right) for designed and post-manufactured FE meshes.

[22,31–33]. These are usually related to the system limitations and a great effort has been made to search for the most adequate processing parameters, as the laser power, the scan speed, or the powder bed preheat temperature [31,33], to reduce the impact of the manufacturing process on the final structure. In this study, the parameters used were optimized for regular orthogonal strut geometries with square or round x-sections originally. These may need to be reviewed for structures with finer details. Another observation is that the differences observed on struts and throat sizes and shape are direction dependent, which indicates that the manufacturing method also induced some anisotropy. This anisotropy induction by the SLS system has also been previously reported [30,31], and it has also been associated with other 3D printing techniques [34]. To analyze how the differences observed between designed and fabricated geometries may influence the final properties (mechanical and permeability), FE meshes of the unit cells were created based on the micro-CT images (here called post-manufacturing meshes), thus presenting geometries closer to the actual ones (Fig. 12), and their computational properties were calculated. When comparing Fig. 12 to Fig. 9, it is possible to notice that these FE meshes are a good representation of the actual unit cells. The permeability values obtained are slightly higher than the designed values (see Fig. 13) and on the other hand, the mechanical properties associated to the post-manufacturing meshes are lower than the prescribed values. These differences on the properties are consistent with higher actual porosity values and smaller struts as well as bigger throat areas on z-direction. In summary, this computational approach can design scaffold architectures for a wide range of situations, as it is possible to simulate not only different strain fields but also to impose minimum

values on the mechanical properties and volume fraction. However, it is important to underline that the manufacturing process used in this study induced some changes on the actual shape and size of the features of the optimized structures and thus their actual properties may not completely match the prescribed values. Overall, this methodology has shown high potential to produce scaffolds for bone tissue engineering and it may originate custom scaffolds for different clinical applications. Funding Portuguese Foundation for Science and Technology (FCT). Competing interest There is not any financial and personal relationship with other people or organizations that could inappropriately influence this work. Acknowledgements Work supported by FCT through the PhD scholarship SFRH/BD/46575/2008 (M. Dias) and Project PTDC/EMEPME/104498/2008. Authors would also like to Krister Svanberg for making the MMA algorithm available for us to use. Appendix A Table A1 Table A1 Designed, STL and manufactured scaffold features.

Porosity (%)

Designed STL Manufactured

Strut size

Designed size (mm) Designed area (mm2 ) STL size (mm) Actual in z-direction

Actual in x/y-direction

Throat size

Designed size (mm) Designed area (mm2 ) STL size (mm) Actual in z-direction

Actual in x/y-direction

2p0

3p0

4p0

54.10 54.46 61.61 ± 1.65

54.10 54.67 56.68 ± 1.15

54.10 54.46 57.94 ± 1.72

Horizontal (mm) Vertical (mm) Area (mm2 ) Horizontal (mm) Vertical (mm) Area (mm2 )

0.80 0.44 0.78 0.56 ± 0.02 0.54 ± 0.03 0.27 ± 0.02 0.59 ± 0.03 0.95 ± 0.02 0.50 ± 0.01

1.20 0.99 1.14 0.96 ± 0.01 0.97 ± 0.005 0.703 ± 0.06 1.22 ± 0.06 1.22 ± 0.08 1.14 ± 0.08

1.60 1.76 1.55 1.41 ± 0.06 1.36 ± 0.03 1.38 ± 0.09 1.43 ± 0.07 1.58 ± 0.05 1.76 ± 0.16

Horizontal (mm) Vertical (mm) Area (mm2 ) Horizontal (mm) Vertical (mm) Area (mm2 )

1.20 1.00 1.22 1.36 ± 0.07 1.40 ± 0.04 1.55 ± 0.16 1.42 ± 0.03 0.99 ± 0.01 1.16 ± 0.02

1.80 2.25 1.86 1.79 ± 0.12 1.86 ± 0.09 2.34 ± 0.10 2.05 ± 0.02 1.74 ± 0.06 2.48 ± 0.04

2.40 4.00 2.45 2.59 ± 0.09 2.60 ± 0.07 4.80 ± 0.30 2.64 ± 0.07 2.30 ± 0.06 4.43 ± 0.14

M.R. Dias et al. / Medical Engineering & Physics 36 (2014) 448–457

References [1] Hutmacher D. Scaffolds in tissue engineering bone and cartilage. Biomaterials 2000;21:2529–43. [2] Yang S, Leong K, Du Z, Chua C. The design of scaffolds for use in tissue engineering: Part I. Traditional factors. Tissue Eng 2001;7(6):679–89. [3] Sanz-Herrera J, Kasper C, van Griensven M, Garcia-Aznar J, Ochoa I, Doblare M. Mechanical and flow characterization of Sponceram carriers: evaluation by homogenization theory and experimental validation. J Biomed Mater Res B: Appl Biomater 2008;87(1):42–8. [4] Babis G, Soucacos P. Bone scaffolds: the role of mechanical stability and instrumentation. Injury 2005;36:38–44. [5] Jones A, Arns C, Hutmacher D, Milthorpe B, Sheppard A, Knackstedt M. The correlation of pore morphology, interconnectivity and physical properties of 3D ceramic scaffolds with bone ingrowth. Biomaterials 2009;30:1440–51. [6] Hollister SJ. Porous scaffold design for tissue engineering. Nat Mater 2005;4:518–90. [7] Mitsak A, Kemppainen J, Harris M, Hollister S. Effect of polycaprolactone scaffold permeability on bone regeneration in vivo. Tissue Eng A 2011;17:1831–9. [8] Sanz-Herrera J, Garcia-Aznar J, Doblare M. A mathematical approach for tissue regeneration inside a specific type of scaffold. Biomech Model Mechanobiol 2007;7:355–66. [9] Jeong C, Zhang H, Hollister S. Three-dimensional poly(1,8-octanediol-cocitrate) scaffold pore shape and permeability effects on sub-cutaneous in vivo chondrogenesis using primary chondrocytes. Acta Biomater 2011;7:505–14. [10] Kemppainen J, Hollister S. Differential effects of designed scaffold permeability on chondrogenesis by chondroyctes and bone marrow stromal cells. Biomaterials 2010;31:279–87. [11] Ochoa I, Sanz-Herrera J, Aznar J, Doblare M, Yunos D, Boccaccini A. Permeability evaluation of 45S5 Bioglass® -based scaffolds for bone tissue engineering. J Biomech 2009;42:257–60. [12] Li S, Wijn J, Li J, Layrolle P, Groot K. Macroporous biphasic calcium phosphate scaffold with high permeability/porosity ratio. Tissue Eng 2003;9:535–48. [13] Chor M, Li W. A permeability measurement system for tissue engineering scaffolds. Meas Sci Technol 2007;18:208–16. [14] Dias M, Fernandes PR, Guedes JM, Hollister SJ. Permeability analysis of scaffolds for bone tissue engineering. J Biomech 2012;45:938–44. [15] Truscello S, Kerckhofs G, Van Bael S, Pyka G, Schrooten J, Van Oosterwyck H. Prediction of permeability of regular scaffolds for skeletal tissue engineering: a combined computational and experimental study. Acta Biomater 2012;8:1648–58. [16] Lee KW, Wang S, Lu L, Jabbari R, Currier B, Yaszemski M. Fabrication and characterization of poly(propylene fumarate) scaffolds with controlled pore structures using 3-dimensional printing and injection molding. Tissue Eng 2006;12(10):2801–11. [17] Bendsøe MP, Sigmund O. Topology optimization theory – methods and applications. Berlin/Heidelberg/New York: Springer; 2003.

457

[18] Kang HS, Long JP, Goldner GDU, Goldstein SA, Hollister SJ. A paradigm for the development and evaluation of novel implant topologies for bone fixation: implant design and fabrication. J Biomech 2012;45:2241–7. [19] Guest J, Prévost J. Design of maximum permeability material structures. Comput Methods Appl Mech Eng 2007;196:1006–17. [20] Sturm S, Zhou S, Mai Y, Li Q. On stiffness of scaffolds for bone tissue engineering – a numerical study. J Biomech 2010;4:31738–44. [21] Zhou S, Li Q. A variational level set method for the topology optimization of steady-state Navier–Stokes flow. J Comput Phys 2008;227:10178–95. [22] Williams JM, Adewunmi A, Schek RM, Flanagan CL, Krebsbach PH, Feinberg SE, et al. Bone tissue engineering using polycaprolactone scaffolds fabricated via selective laser sintering. Biomaterials 2007;2:4817–27. [23] Coelho PG, Fernandes PR, Rodrigues HC, Cardoso JB, Guedes JM. Numerical modeling of bone tissue adaptation – a hierarchical approach for bone apparent density and trabecular structure. J Biomech 2009;42(7):830–7. [24] Guedes J, Kikuchi N. Preprocessing and post processing for materials based on the homogenization method with adaptive finite element method. Comput Methods Appl Mech Eng 1990;83:143–98. [25] Becker EB, Carey GF, Oden JT. Finite elements: an introduction, vol. I. New Jersey: Prentice-Hall; 1981. [26] Hollister SJ, Lin CY. Computational design of tissue engineering scaffolds. Comput Methods Appl Mech Eng 2007;196:31–2. [27] Terada K, Ito T, Kikuchi N. Characterization of the mechanical behaviors of solid–fluid by the homogenization method. Comput Methods Appl Mech Eng 1998;153:223–57. [28] Svanberg K. The method of moving asymptotes – a new method for structural optimization. Int J Numer Methods Eng 1987;24:359–73. [29] Coelho PG, Fernandes PR, Rodrigues HC. Multiscale modeling of bone tissue with surface and permeability control. J Biomech 2011;44(2):321–9. [30] Cahill S, Lohfeld S, McHugh PE. Finite element predictions compared to experimental results for the effective modulus of bone tissue engineering scaffolds fabricated by selective laser sintering. J Mater Sci Mater Med 2009;20(6):1255–62. [31] Eshraghi S, Das S. Mechanical and microstructural properties of polycaprolactone scaffolds with 1-D, 2-D, and 3-D orthogonally oriented porous architectures produced by selective laser sintering. Acta Biomater 2010;6(7):2467–76. [32] Eosoly S, Brabazon D, Lohfeld S, Looney L. Selective laser sintering of hydroxyapatite/poly-e-caprolactone scaffolds. Acta Biomater 2010;6: 2511–7. [33] Partee B, Hollister SJ, Das S. Selective laser sintering process optimization for layered manufacturing of CAPA® 6501 polycaprolactone bone tissue engineering scaffolds. J Manuf Sci Eng: Trans ASME 2006;128:531–40. [34] Butscher A, Bohner M, Hofmann S, Gauckler L, Müller R. Structural and material approaches to bone tissue engineering in powder-based three-dimensional printing. Acta Biomater 2011;7(3):907–20.