2.1 Bone ingrowth simulation
The simulation was used to examine bone growth during the rehabilitation phase, which was numerically implemented using coding the subroutine VUMAT of the commercial finite element software Abaqus (Dassault Systèmes, USA).
2.1.1 Material
The choice of the scaffold material is a crucial aspect in the design of a bone implant. Although biodegradable scaffolds have been shown to stimulate tissue growth and elicit minimal immune response after they dissolve, permanent biocompatible metals remain the most commonly used materials due to their strong mechanical strength, ease of large-scale manufacturing, and relatively lower cost. For this study, Ti-6Al-4V was selected as the scaffold material for all architectural designs. This material exhibits a Young's modulus of 107 GPa, a Poisson's ratio of 0.3, and a density of 4500 kg/m3. The scaffold pores are initially filled with interstitial fluid (ISF), which will be partially replaced by bone tissue over time. The materials in the model, including the scaffold, bone, and ISF, were assumed to be isotropic and linearly elastic. Additionally, the ISF was assumed to be nearly incompressible.
2.1.2 Scaffold model design
The bone scaffold model was designed in Solidworks® and imported into Abaqus for further simulation. Each scaffold unit cell is a cube that is 1 millimetre on each side. The cubic and circular designs, as shown in Fig. 1a and b respectively. Figure 1c indicates FD-cube (face-diagonal cube) unit cell, and Fig. 1d shows V-Octet (Void Octet) design. The novel-designed scaffold structure Hexnaoid is manifested in Fig. 1e.
2.1.3 Meshing and voxelization
The computational domain was Hex meshed into identical cubic elements with a side length of 50 µm. The titanium scaffold and the porous cavity were allocated as two separate sets. Both sets were converted from the continuous domain into a discrete grid block through outsourced voxelization scripts. The voxelization was based on the scaffold’s structure, as indicates in Fig. 1f-j. Each element can have one of the following two states for the scaffold cavity: bone or ISF. In order to display the states during the whole process, we defined a characteristic function χ: If χ = 1, the element was in the state of the scaffold; if χ = 2, the element was in the state of bone; and the element was in the state of ISF when χ = 3.
2.1.4 Boundary conditions
The simulation setup was designed to imitate the loads that would be applied to a scaffold implant in a real-life scenario. The top surface of the scaffold model was subjected to a uniformly distributed pressure, simulating the loads on the model. To ensure homogeneous deformation of the entire model in the z-direction, the bottom surface of the scaffold was fixed and a rigid plate was attached to the top surface. The applied pressure was defined as a trapezoidal pulse with a period of 1 day (see Fig. 2), including relaxation, rising, hold, and fall phases. The latter three phases define the duration of movement, with the relaxation phase denoting idleness. The rising and falling phases of the loading history were set to 0.05 days to avoid sudden changes in the loading history, which could lead to erroneous simulations. The simulation was executed over a duration of 100 days, which is equivalent to the average recovery time of an individual following a fracture (approx. 16 weeks) (Nikolaou, Efstathopoulos, Kontakis, Kanakaris, & Giannoudis, 2009). It was used to examine bone growth during the rehabilitation phase. According to Shefelbine, where the compressive stress applied to the surface of the rigid upper plate was 3 MPa (Shefelbine, Augat, Claes, & Simon, 2005).
2.1.5 Bone remodelling theory
The simulation assumed that bone formation and resorption are primarily facilitated by osteoblasts and osteoclasts, respectively, at the cellular level on the bone surface. The osteogenesis theory posits that these cells can only respond to mechanical signals when they are attached to a scaffold or on the surface of previously formed bone that has undergone differentiation from cartilage(Horwitz & Parsons, 1999). As a result, bone resorption and formation are not considered to occur within the interstitial fluid. To simplify the process, the simulation relied on the locally heterogeneous Strain Energy Density (SED) ψ as the sole stimulus affecting osteogenesis and resorption factors. Based on Husikes theory and Schulte's work (Huiskes et al., 1987; Schulte et al., 2013) the bone remodelling rate u(ψ) represents the thickness change of formed/absorbed bone per unit time, as shown in Fig. 3, and is expressed as:
In the formula, c is a constant and represents how fast bone formation and resorption rates reach the maximum growth rate umax. And ψupper and ψlower are bone formation and resorption thresholds, respectively. The “lazy zone” presents the specific mechanical stress interval in which osteoblastic and osteoclastic cellular activities reach dynamic balance. Between this interval, there is no bone formation or degeneration macroscopically. For each random element n, its local SED ψ (xn) is affected by its neighbouring element m within a sensitive distance D. Within this range, closer element mth has a greater influence on element nth, and vice versa. The formulas of local strain energy density can be formulated as: (Schulte et al., 2013)
where q is the number of other elements that has influences to this element. SED (xm) is the strain energy density of the mth element, and d (xm −xn) is the distance between element n and m. The rate of bone volume fraction αb(t) of a bone element (equal to relative density \({\stackrel{-}{\rho }}_{b}\)(t)3) in the dynamic process depends on the bone remodelling rate u(ψ), and can be defined as the following formula:
where l is the side length of the elements. The Young’s modulus of bone element at time t, Eb(t), was computed by a modified modulus-density relationship, where Eb and EISF are the Young’s moduli of the natural bone and ISF, respectively.
Each element being ossified in a certain degree of bone maturation can contribute to the mechanical properties of the scaffold-bone system only when the relative density \({\stackrel{-}{\rho }}_{b}\)(t)3 is above a certain threshold (Schulte et al., 2013). In other words, when the degree of ossification is low, this element is still considered ISF. Therefore, the judgement criteria of ossification is element relative density \({\stackrel{-}{\rho }}_{b}\)(t)3. When αb(t) of an element is less than a threshold ρthre, the element is changed into ISF (resorption), i.e., χ from 2 to 3. On the contrary, when density \({\stackrel{-}{\rho }}_{b}\)(t)3 of an element is greater than ρthre, the element undergoes osteogenesis (formation), i.e., χ from 3 to 2.
2.1.5 Input parameters
In this instance, titanium made up the scaffold; its mechanical characteristics are already mentioned in section 2.1.2. Frost measured that the maximum rate of bone resorption or production during bone remodelling was 2 mm3/mm2/year, which translated into umax=0.005 mm/day in the current simulations (Frost, 1990). The lower and upper criteria (ψlower and ψupper) were modified based on literature (Schulte et al., 2013). As indicated in Section 2.1, the boundary conditions for the rehabilitative exercise level were 3 MPa. Table 1 contains a list of all input parameters utilised in the simulation.
Table 1
Summarization of the bone remodelling theory parameters used in the FEA simulation
Parameters
|
|
Value
|
Unit
|
Constant
|
c
|
0.5 (Schulte et al., 2013)
|
Mm/Mpa/day
|
Maximal formation/resorption velocity
|
umax
|
0.005 (Frost, 1990)
|
Mm/day
|
Resorption threshold
|
Ψlower
|
0.01 (Schulte et al., 2013)
|
MPa
|
Formation threshold
|
Ψupper
|
0.02 (Schulte et al., 2013)
|
MPa
|
Influence distance
|
D
|
52 (Schulte et al., 2013)
|
µm
|
Young’s modulus of mature bone
|
Eb
|
20 (Morgan, Bostrom, & van der Meulen, 2015)
|
GPa
|
Poisson’s modulus of scaffold and bone
|
ν
|
0.3
|
-
|
Young’s modulus of ISF
|
EISF
|
0.01
|
GPa
|
Poisson’s modulus of ISF
|
νISF
|
0.49
|
-
|
State change threshold
|
ρthre
|
0.01
|
--
|
2.2 Scaffold mechanical compression test and simulation
The mechanical strength of the Hexanoid scaffold was evaluated through a combination of experimental quasi-static compression testing and numerical simulation using finite element analysis (FEA) software. The architecture and cross-section of the Hexanoid scaffold used in the mechanical compression tests are depicted in Figure 4.a. The experimental mechanical compression test was conducted three times using 3D-printed specimens composed of a 4×4×6 unit cell structure, made of Ti-6Al-4V alloy (as depicted in Figure 4.b.). On the other hand, the module used in the FEA simulation was designed as a 3×3×3 unit cell structure to reduce simulation time and complexity while ensuring the minimization of simulation error, as shown in Figure 4.c. (Simsek, Akbulut, Gayir, Basaran, & Sendur, 2021). Regarding the boundary condition, the lower surface of the compression cylinder was constrained in all direction, and a displacement of 0.15mm (5% of strain) was applied to the upper cylinder surface. The corresponding reaction force to the displacement was recorded as output. A stress vs. strain curved was plotted considering the loaded surface area, and Young’s modulus of the elastic region was calculated. The same simulation process and principle were applied in all other scaffold designs to compare their mechanical strength.