Computational Analysis of Bone Remodeling in the Proximal Tibia Under Electrical Stimulation Considering the Piezoelectric Properties

The piezoelectricity of bone is known to play a crucial role in bone adaptation and remodeling. The application of an external stimulus such as mechanical strain or electric field has the potential to enhance bone formation and implant osseointegration. Therefore, in the present study, the objective is to investigate bone remodeling under electromechanical stimulation as a step towards establishing therapeutic strategies. For the first time, piezoelectric bone remodeling in the human proximal tibia under electro-mechanical loads was analyzed using the finite element method in an open-source framework. The predicted bone density distributions were qualitatively and quantitatively assessed by comparing with the computed tomography (CT) scan and the bone mineral density (BMD) calculated from the CT, respectively. The effect of model parameters such as uniform initial bone density and reference stimulus on the final density distribution was investigated. Results of the parametric study showed that for different values of initial bone density the model predicted similar but not identical final density distribution. It was also shown that higher reference stimulus value yielded lower average bone density at the final time. The present study demonstrates an increase in bone density as a result of electrical stimulation. Thus, to minimize bone loss, for example, due to physical impairment or osteoporosis, mechanical loads during daily physical activities could be partially replaced by therapeutic electrical stimulation.

Introduction

The human skeletal system consists of bones and joints, which maintains the structural integrity of the body, provides sites for muscle attachment, and facilitates body movements (Cowin, 2001). Bone adapts its structure in response to changes in its mechanical loading environment and this adaptation process is known as bone remodeling (Robling and Turner, 2009). This process has a significant impact on the individual’s health and thus bone remodeling study is of prime importance. Additionally, piezoelectricity plays a vital role in bone adaptation and remodeling processes (Mohammadkhah et al., 2019). Therefore, in order to achieve better understanding of mechanical and electrical interactions that occur during these processes, computational analysis of piezoelectric bone remodeling is of great interest in musculoskeletal biomechanics.

In the literature, many mathematical models of bone remodeling are based on the qualitative observations of Wolff (1893), and these models have been implemented using the finite element method to simulate the bone response to mechanical loading (Weinans et al., 1992; Fernández et al., 2010; Fernández et al., 2012a; Garzón-Alvarado et al., 2012; Mauck et al., 2016; Garijo et al., 2017; Quilez et al., 2017) and have been reviewed in (Gerhard et al., 2009). The bone remodeling algorithms proposed in the past are based on strain (Wang et al., 2016), stress (Gong et al., 2013), strain energy density (SED) (Weinans et al., 1992; Fernández et al., 2010; Bouguecha et al., 2011), deformation (Papathanasopoulou et al., 2002), and mechanical damage (Garcia-Aznar et al., 2005; Hambli et al., 2011) or a combination of these, which elucidate many aspects of bone adaptation. The rate of change in bone density can be linear or nonlinear with respect to the applied mechanical load. Moreover, the increment and decrement in this rate can be similar or different (Su et al., 2019). Although bone remodeling has multiple aspects, these are confined to study the response of bone to a particular loading type.

Piezoelectricity can explain the involvement of electrical signals and mechanical loads in the bone adaptation process. Yasuda (1954) was the first one to report the direct and converse piezoelectric properties of the dry bone and these findings have been confirmed by other investigators (Fukada and Yasuda, 1957; Johnson et al., 1980). They attributed the piezoelectric behavior of bone to matrix piezoelectricity, in which application of mechanical shear to collagen fibers generates electrical charges. It was found that the quasi-hexagonal symmetry of the collagen fiber is responsible for the shear piezoelectric effect in bone (Minary-Jolandan and Yu, 2009). Based on experimental observations (Pienkowski and Pollack, 1983; Hastings and Mahmud, 1988), the streaming potential causes piezoelectricity in wet bone and is the main source for strain-generated potentials (SGPs). Streaming potential arises from the fluid flow through charged surfaces. Mechanical deformation causes fluids to flow through the canaliculi and correspondingly ions flow is generated against the oppositely charged walls, leading to a potential difference between two points along the stream (Mohammadkhah et al., 2019). Until now, the specific mechanisms for bone piezoelectricity remain unclear. However, there are few computational models of bone remodeling that take the piezoelectricity of bone into account (Qu et al., 2006; Fernández et al., 2012a; Garzón-Alvarado et al., 2012; Beheshtiha and Nackenhorst, 2015; Cerrolaza et al., 2017; Duarte et al., 2018) and have been summarized in recent review by Mohammadkhah et al. (2019).

Several studies have described the piezoelectric effects on bone biophysics. The linear piezoelectric theory for bone remodeling was presented by Gjelsvik (1973). It was observed in experimental studies that bone possesses electric and dielectric properties, which are frequency-dependent (Bassett and Becker, 1962; Bassett et al., 1964; Su et al., 2016) and these properties are important for the hypothesized feedback mechanism for the bone remodeling process and the therapeutic electrical stimulation for bone healing process (Ramtani, 2008). It has been demonstrated that the electromagnetic field affects the bone remodeling and healing process under the influence of mechanical and electrical loadings (Qu et al., 2006; Qu and Yu, 2011). Also, there are models that considered bone remodeling under the combined action of mechanical, electrical, and thermal loads (Qin, 2003; Qin et al., 2005). However, the effect of these stimuli on the remodeling process is still ambiguous. For example, an understanding of this could be useful for improving the implant design (Lian et al., 2014). Although there are many femoral bone remodeling studies, similar studies on tibial bone are very few (Perez et al., 2010; Robalo, 2011; Fang et al., 2013; González-Carbonell et al., 2015). Moreover, there is no investigation in the literature on bone remodeling in the human tibia considering its piezoelectric properties. Therefore, the aim of our present study is to investigate piezoelectric bone remodeling in the human tibia under electrical stimulation. This study contributes to a better understanding of tibial bone response to electromechanical loads that can aid in designing a protocol for therapeutic electrical stimulations to minimize bone loss, for example, due to physical impairment or osteoporosis.

Material and Methods

Mathematical Formulation of the Piezoelectric Bone Remodeling

Several bone remodeling algorithms have been proposed to compute the evolution of bone density under mechanical loadings; however, only a few of these algorithms have been considering bone piezoelectricity. Based on the piezoelectric strain-adaptive bone remodeling algorithm proposed by Fernández et al. (2012a), let $\Omega$ be an open-bounded domain (see Figure 1A), and $\Gamma$=$\partial \Omega$ be its boundary. This boundary was considered to be Lipschitz continuous. It has been split into a Neumann boundary${\Gamma }_N$ and a Dirichlet boundary ${\Gamma }_D$. The density of volume forces acting in domain $\Omega$ is denoted by ${\mathit{f}}_{\mathit{B}}$, and the density of traction forces ${\mathit{f}}_{\mathit{N}}$ was applied on${\Gamma }_N$. It has been assumed that the bone is fixed at ${\Gamma }_D$, i.e., here $\mathit{u}$ = 0 holds for the displacement vector. Next, the density of volume electric charges present in domain$\Omega$ is represented by $q_B$ and the density of surface electric charges applied externally on${\Gamma }_N$ is represented by $q_N$. Let $\phi$ be the electric potential. An electric potential ${\phi }_D=0$ was applied to the fixed boundary (Fernández et al., 2012a). Let [0,$T$], $T$ > 0 be the period of interest and $\mathit{\nu }\left(\mathit{x}\right)$ be the outward unit normal vector to $\Gamma$ at a point $\mathit{x}$. Here, the same decomposition of the boundary has been used for imposing the boundary conditions for the mechanical displacement and the electric potential (Fernández et al., 2010). Note that the symbols in bold represent vectors, tensors, or matrices.

FIGURE 1

FIGURE 1. (A) The tibia bone domain and schematic representation of boundary conditions implemented for piezoelectric bone remodeling simulations under electrical stimulation. (B) Exemplary representation of the direct and converse piezoelectric effect in tibia bone.

The linearized strain tensor $\mathit{\epsilon }\left(\mathit{u}\right)$ can be defined as

${\epsilon }_{ij}\left(\mathit{u}\right)=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\text{,}i,j=1\text{,}\dots ,d,(1)$

where $\mathit{u}$ represents the displacement field and d is the order of symmetric matrices (3×3). Similar to Weinans et al. (1992), the bone was assumed to be isotropic and linear elastic. For the stress $\mathit{\sigma }$, the constitutive law can be written as follows:

$\mathit{\sigma }=\mathit{\sigma }\left(\mathit{u}\right)=2\mu \left(\rho \right)\mathit{\epsilon }\left(\mathit{u}\right)+\lambda \left(\rho \right)\text{Div}\left(\mathit{u}\right)ℐ-\alpha \left(\rho \right){\mathit{ℰ}}^{\ast }\mathbf{E}\left(\phi \right)\text{in}\overline{\mathit{\Omega }}\times \left[0,\mathit{T}\right],(2)$

where $\mu \left(\rho \right)$ and $\lambda \left(\rho \right)$ are Lamé coefficients and were considered to be dependent on the bone apparent density $\rho$, $\text{Div}$ denotes the divergence operator, and ℐ represents the identity operator (Fernández et al., 2010). For the plane strain condition or the three-dimensional model, Lamé coefficients can be defined in terms of elastic modulus $E\left(\rho \right)$and Poisson’s ratio $k\left(\rho \right)$ as follows:

$\mu \left(\rho \right)= \frac{E\left(\rho \right)}{2\left(1+k\left(\rho \right)\right)}\mathrm{a}\mathrm{n}\mathrm{d}\lambda \left(\rho \right)= \frac{k\left(\rho \right)E\left(\rho \right)}{1-k^2\left(\rho \right)}.(3)$

The Poisson’s ratio was considered to be not dependent on ρ (and thus, k(ρ) = k). The following equation was employed for elastic modulus depending on the apparent bone density:

$E\left(\rho \right)=M{\rho }^{\gamma },(4)$

where $M$ and $\gamma$ are positive constitutive constants that characterize the behavior of bone (Weinans et al., 1992). Similar to the elastic modulus, a constitutive function$\alpha \left(\rho \right)$ was assumed to be dependent on the apparent bone density function and expressed as (Fernández et al., 2012a),

$\alpha \left(\rho \right)= {\rho }^{\gamma }.(5)$

Further, ℰ^∗ denotes the transpose of the third-order piezoelectric tensor ℰ described below. As a conservative field, the stationary electric field E can be computed from the gradient of the electric potential $\phi$ (van Rienen, 2001):

$\mathit{E}=-\nabla \phi .(6)$

The following first-order ordinary differential equation was employed to calculate the evolution of the apparent bone density function (Weinans et al., 1992),

$\frac{d\rho }{dt}=B\left(\frac{U\mathit{ }\left(\mathit{\sigma }\left(\mathit{u}\right), \mathit{\epsilon }\left(\mathit{u}\right)\right)}{\rho }- S_r\right)\text{in} \mathit{\Omega }\times \left(0,T\right),(7)$

where the values of experimental constants $B$ and $S_r$ are mentioned in Table 1. Further, the strain energy density (SED) as mechanical stimulus $U(\mathit{\sigma }\left(\mathit{u}\right), \mathit{\epsilon }\left(\mathit{u}\right))$ can be given as:

$U\left(\mathit{\sigma }\left(\mathit{u}\right),\mathit{\epsilon }\left(\mathit{u}\right)\right)=\frac{1}{2}\mathit{\sigma }\left(\mathit{u}\right) :\mathit{\epsilon }\left(\mathit{u}\right),(8)$

where “$:$” denotes the inner product and the apparent density has been considered to be restricted as follows:

${\rho }_a \le \rho \le {\rho }_b,(9)$

where the minimum density ${\rho }_a$ and the maximum density ${\rho }_b$ correspond to the resorbed bone and the cortical bone, respectively. The constitutive law for the electric displacement $\mathbf{D}$ can be expressed as:

$\mathbf{D}= \text{α}\left(\rho \right)\mathit{ℰ}\mathit{\epsilon }\left(\mathit{u}\right)+\alpha \left(\rho \right)\mathit{\beta }\mathbf{E}\left(\phi \right),(10)$

where $\mathit{\beta }$ is the electric permittivity tensor (Batra and Yang, 1995). The constitutive equations for the stress $\mathit{\sigma }$ (Eq. 2) and the electric displacement $\mathbf{D}$ (Eq. 10) govern the piezoelectric behavior of bone. When subjected to mechanical loading, bone generates an electric charge (direct piezoelectric effect), and conversely, when an electrical charge is applied, strains/stresses can appear in bone (converse piezoelectric effect) (Fernández et al., 2012a; Fernández et al., 2012c) (see Figure 1B). Similar to other studies (Fotiadis et al., 1999; Qin and Ye, 2004; Fernández et al., 2012a, Fernández et al., 2012c; Duarte et al., 2018), the bone was assumed to behave like a crystal with hexagonal symmetry, i.e., the third-order piezoelectric stress tensor ℰ is defined by four values and the electric permittivity tensor (dielectric tensor) $\mathit{\beta }$ is a diagonal matrix with two constants. These tensors can be written in the following matrix form:

$\mathit{ }\mathit{ℰ}\mathit{=}\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ e_{31}& e_{31}& e_{33}\end{array} \begin{array}{ccc}{e}_{14}& e_{15}& 0\\ e_{15}& -e_{14}& 0\\ 0& 0& 0\end{array}\right)\text{and}\mathit{\beta }=\left(\begin{array}{ccc}{\beta }_{11}& 0& 0\\ 0& {\beta }_{11}& 0\\ 0& 0& {\beta }_{33}\end{array}\right),(11)$

where the third principal direction represents the longitudinal direction of the tibia bone (Fernández et al., 2012c). Forces resulting from the musculoskeletal system (e.g., muscle forces) were applied as Neumann boundary conditions to the medial and lateral condyles. Further, the following coupled linear variational equations were solved to calculate the displacement field $\mathit{u}$ and the electric potential $\phi$.

${\displaystyle \underset{\Omega }{\int }2\mu \left(\rho \right)\mathit{\epsilon }\left(\mathit{u}\right):\mathit{\epsilon }\left(\mathit{v}\right)+ \lambda \left(\rho \right)\text{Tr}\left(\mathit{\epsilon }\left(\mathit{u}\right)\right)\text{Tr}\left(\mathit{\epsilon }\left(\mathit{v}\right)\right)d\mathit{x}}={\displaystyle \underset{\Omega }{\int }{\mathit{f}}_B\left(t\right)\cdot \mathit{v}d\mathit{x}}$

$+ {\displaystyle \underset{{\Gamma }_N}{\int }{\mathit{f}}_N\left(t\right)\cdot \mathit{v}d\Gamma }-{\displaystyle \underset{\Omega }{\int }\left(\rho {\left(t\right)}^{\gamma }{\mathit{ℰ}}^{\ast }\nabla \phi \left(t\right),\mathit{\epsilon }\left(\mathit{v}\right) \right)d\mathit{x}},(12)$

${\displaystyle \underset{\Omega }{\int }\left(\rho {\left(t\right)}^{\gamma } \mathit{\beta } \nabla \phi \left(t\right), \nabla \psi \right)d\mathit{x}}={\displaystyle \underset{\Omega }{\int }{q}_B\left(t\right) \psi d\mathit{x}}$

$+ {\displaystyle \underset{{\text{Γ}}_N}{\int }{q}_N\left(t\right) \psi d\Gamma }+ {\displaystyle \underset{\Omega }{\int }\left(\rho {\left(t\right)}^{\gamma }\mathit{ℰ}\left(\mathit{u}\left(t\right)\right),\nabla \psi \right)d\mathit{x}},(13)$

where $\text{Tr}$ denotes the classical trace operator, $\mathit{v}$ and $\psi$ are the test functions, $d\mathit{x}$ denotes the differential element for integration over domain $\mathit{\Omega }$, and $q_B=\text{div }\mathbf{D}$. For simplicity, inertial effects were neglected and the bone remodeling formulation was restricted to isothermal and quasi-static conditions. More details about this model can be found in (Fernández, 2010; Fernández et al., 2010; Fernández et al., 2012a).

TABLE 1

TABLE 1. Values of the parameters used in the finite element analyses.

For the numerical implementation of this formulation in the Python-based open-source finite element software FEniCS (www.fenicsproject.org, version 2019.1.0, GNU-GPL) (Logg et al., 2012; Alnæs et al., 2015), a mixed-function space (vector function space for the mechanical displacement and scalar function space for the electric potential) was used with Lagrange elements of order 2. Throughout the study, the time derivatives were discretized using the explicit fourth-order Runge-Kutta method.

${\rho }_{n+1}={\rho }_n+\frac{\Delta t}{6} \left(k_1 +2k_2 +2k_3 + k_4\right),(14)$

$k_1 =f\left(t_n, {\rho }_n\right),(15)$

$k_2 =f\left(t_n+\frac{\Delta t}{2}, {\rho }_n+\Delta t\frac{k_1}{2}\right),(16)$

$k_3 =f\left(t_n+\frac{\Delta t}{2}, {\rho }_n+\Delta t\frac{k_2}{2}\right),(17)$

$k_4 =f\left(t_n+\frac{\Delta t}{2}, {\rho }_n+h{k}_3\right),(18)$

where Δt is the size of the time increment, ${\rho }_{n+1}$ and ${\rho }_n$ represent the bone density for the new and the current time step, respectively, and $f\left(t_n,{\rho }_n\right)=B\left(\frac{\mathit{U}(\mathit{\sigma }\left(\mathit{u}\right),\mathit{\epsilon }\left(\mathit{u}\right))}{{\rho }_n}-S_r\right)$. The values of the parameters used in the following simulations are listed in Table 1.

Setting up a Finite Element Model Using the Open-Source Framework

In this study, open-source software were used to ensure reproducibility of the simulation results, which is important for the evaluation of scientific work. Within our implementation of the open-source framework (Figure 2A), the Python-based open-source software packages used were ITK-SNAP (http://www.itksnap.org/, version 3.6.0, GNU-GPL) (Yushkevich et al., 2006) for segmentation, Salome (https://www.salome-platform.org/, version 8.5.0, GNU-LGPL) (Ribes and Caremoli, 2007) for designing and meshing, Gmsh (http://gmsh.info/, version 4.4.1, GNU-GPL) (Geuzaine and Remacle, 2009) and command-line tool dolfin-convert (GNU-LGPL) for pre-processing, the finite element software FEniCS for solving and post-processing, and Paraview (https://www.paraview.org/, version 5.0.1, 3-Clause BSD License) (Ahrens et al., 2005) for visualizing the simulation results. More details of the steps followed to set up the finite element model are published in (Bansod et al., 2021; Bansod and van Rienen, 2019). Figure 2B outlines the implemented bone remodeling algorithm, where the material properties were updated in each iteration through the positive feedback loop. The piezoelectricity of bone was incorporated into the algorithm by considering the matrix piezoelectricity (see Eq. 11) (Fernández et al., 2012a; Mohammadkhah et al., 2019). Consequently, the bone density distribution predicted by the finite element model is for dry bone. Following Open Science principles, the source code developed in this study is available at https://github.com/YDBansod/Bone_Remodelling.

FIGURE 2

FIGURE 2. (A) Open-source software framework used to simulate tibial bone remodeling. (B) Flowchart of the piezoelectric bone remodeling algorithm implemented in the open-source finite element software FEniCS.

Loads and Boundary Conditions

This study focuses on tibial bone remodeling under electrical stimulation, which has not yet been investigated. The tibia geometry was meshed using 10-node isoparametric tetrahedral elements. The mesh includes 4,132 elements and 1,507 nodes (see Figure 3A). The end-to-end distance between the medial and lateral condyles W is 71.10 mm and the diameter of the diaphysis D is 21.90 mm. The height of the proximal tibia L is 108.68 mm. Here, the lateral condyle-plateau is approximately 1.68 times that of the medial condyle-plateau. An obvious drawback of this two-dimensional model is the lack of connection between the cortical diaphysis regions. As a solution to this problem, an additional side-plate was considered (Figure 3B) joining these regions similar to that of the two-dimensional femur (Weinans et al., 1992; Fischer et al., 1997; Fernández et al., 2010; Garijo et al., 2014; Bansod et al., 2021). The cortical regions were joined only at medial and lateral nodes highlighted with black dots (see Figure 3C). The side-plate P-Q-R-S consists of 1,439 and 562 elements and nodes, respectively, and had mechanical properties analogous to the cortical bone. Additionally, the remodeling capacities of this plate were restricted.

FIGURE 3

FIGURE 3. Finite element meshed models: (A) the human proximal tibia and (B) the side-plate. (C) Both tibia and side-plate were connected to each other at medial nodes (on the edge P-Q) and lateral nodes (on the edge S-R) highlighted by black dots. Note that the electrical boundary conditions are highlighted in orange. (D) Lateral view showing the thickness of the proximal tibia.

Two Dirichlet boundary conditions were imposed at the lower region of the tibia, where the bottom face was restrained in the vertical direction, and to prevent rigid body movement, it was centrally constrained in both vertical and horizontal directions (see Figure 3C) (Nyman et al., 2004). Additionally, this central edge was electrically grounded (electric potential = 0 V). To perform multiplication of tensors of different orders, the tibia and the side-plate were modelled as a slice with a uniform thickness of 1 and 0.1 mm, respectively. Since partial differential equation (PDE) solver FEniCS does not support hexahedral mesh, both domains were meshed using tetrahedral elements with one element along the thickness (see Figure 3D). The anterior and posterior faces of the tibia and the side-plate were constrained along the z-axis to impose plane strain conditions. The total time of simulation and step size were set as t = 300 days (i.e., remodeling period) and Δt = 0.125 days, respectively.

Forces resulting from muscle activities were applied as Neumann boundary conditions to the medial and lateral condyles. Mechanical loading conditions corresponding to the walking activity were simulated through the joint reaction force (approximately 3 times bodyweight of 70 kg) at the condylar region (Duda et al., 2001; Nyman et al., 2004; Perez et al., 2010; Quilez et al., 2017). The walking activity was considered to be represented by three cyclic loading scenarios. In the first loading scenario, the joint reaction force was distributed equally to each condyle in a vertical direction (Figure 4A). In the second loading scenario, the joint reaction force was distributed 70% over the medial condyle and 30% over the lateral condyle, while in the third loading scenario it was distributed 30 and 70% over the medial and lateral condyles, respectively. In the second and third loading scenarios, the joint reaction force was inclined 5^∘ to the vertical direction such that its horizontal-component was directed towards the medial region (see Figures 4B,C) (Nyman et al., 2004; Perez et al., 2010; Quilez et al., 2017). Each loading scenario consists of a pair of parabolic distributed loads acting over the condyles and the load values considered are presented in Table 2. These cyclic loading scenarios with different frequencies were applied independently and consecutively (Figure 4D), where each iteration represents a day. The simulation was initially run for t = 300 days and with the obtained density distribution as an initial configuration, it was further assumed that between 300 and 400 days, the daily physical activity of the person was reduced due to injury or post-surgery. To incorporate this, the mechanical loads were applied once every 4 days between days 300–400 and referred to as “reduced physical activity” by Fernández et al. (2012a). Similar to piezoelectric material, when electrical stimulation is applied to bone, a corresponding mechanical displacement is obtained, which in this case leads to a change in bone density. Accordingly, during the reduced physical activity period (i.e., days 300–400), a surface electric charge q_N = 4 × 10⁻⁹ C/mm² was applied to the lower end of the medial condyle (see Figure 4E).

FIGURE 4

FIGURE 4. Loading and boundary conditions: (A) loading scenario 1 (B) loading scenario 2, and (C) loading scenario 3. Mark that for loading scenarios 2 and 3, the joint reaction force is inclined by 5^∘ to the vertical direction. (D) Loading pattern applied independently and sequentially. (E) Therapeutic electrical stimulation in the form of surface electric charge applied to the lower end of the medial condyle during the period of reduced physical activity. Note that the electrical boundary conditions are highlighted in orange.

TABLE 2

TABLE 2. Values of the forces applied on each condyle (Perez et al., 2010; Quilez et al., 2017).

The bone remodeling simulations started with a uniform initial density and by imposing the relevant boundary conditions (mechanical and electrical), the evolution of bone density was calculated. Here, the uniform initial bone density ${\rho }_0$ of 0.8 g/cm³ was selected because it is the mean of the bone density values (minimum 0.001 g/cm³ and maximum 1.74 g/cm³) used in the simulations (see Table 1). Starting from homogeneous density distribution, bone changes its shape in response to the applied external loads and at the end of the simulation, shows heterogeneous bone density distribution. In the present study, the bone density evolution was computed using the element-based approach (i.e., here the density was assumed to be constant elementwise). Further, the bone was considered to be homogeneous and isotropic linear-elastic material. Small displacement theory was assumed throughout this study (i.e., the displacement is so small that the deformed configuration of bone is not much different from the original shape). A direct LU solver was used to solve the linear asymmetric system. The simulations were performed on a computer with an Intel(R) Xeon(R) CPU E5-2687 W v6 @ 3.10GHz, 256GB RAM, and eight physical cores using FEniCS and it took approximately 12 s for each time-iteration.

Parametric Study

Since the implemented bone remodeling algorithm is complex, iterative, and includes several parameters, a detailed parametric study is necessary. The influence of model parameters such as the uniform initial bone density ${\rho }_0$ and the reference stimulus S_r on the final bone density distribution was investigated. Concerning the initial bone density, five different values within the range of 0.2–1.4 g/cm³ and for the reference stimulus, six values within the range of 0.002–0.008 J/g were considered. Additionally, the variations within a range of ±1%, ±3%, and ±5% of the selected initial bone density of 0.8 g/cm³ were also considered to evaluate the model performance. The other model parameters were kept unchanged while performing these studies. The time step and mesh convergence studies were also conducted, and optimum mesh size of 1 mm and time step of 0.125 days were used in all simulations.

Conversion of Hounsfield Units Into Bone Mineral Density

This study is based on a fresh human tibia specimen from a female donor (58 years old and 163 cm height) without any history of orthopedic injury or surgery. The specimen was scanned using CT (Aquilion 32, Toshiba, Neuss, Germany) and the images were saved in the digital imaging and communications in medicine (DICOM) format. By importing these images into AMIRA^® software (v.5.4.1, Zuse Institute Berlin, Berlin, Germany; Thermo Fisher Scientific, Waltham, MA, United States), the bone surface was segmented. In the literature, there are several relations relating HU to apparent bone density (Wirtz et al., 2000; Peng et al., 2007; Knowles et al., 2016). In the present study, to associate the HU values with the bone density at every point i of the tibia bone, the equation from Garijo et al. (2017) was modified as follows:

${\rho }_i={\rho }_a+ \frac{{\rho }_b-{\rho }_a}{H{U}_{max}- H{U}_{min}} \left(H{U}_i- H{U}_{min}\right).(19)$

where the HU maximum ($H{U}_{max}$) and minimum ($H{U}_{min}$) values were obtained from the CT data and correlated with bone density values of ${\rho }_b$ = 1.74 g/cm³ (cortical bone) and ${\rho }_a$ = 0.01 g/cm³ (resorbed bone), respectively (Weinans et al., 1992; Fernández et al., 2012a). The study was approved by the local ethical committee (registration number A2017-0110). Concerning the validation of simulation results, the predicted bone density distributions were qualitatively and quantitatively compared with the values obtained from the CT data. More precisely, for qualitative analysis, the obtained density distribution was visually matched with the CT image of the same tibia. Additionally, for quantitative analysis, the root mean square (RMS) error and mean deviation (MD) were calculated using the absolute differences between the bone densities estimated from the simulation and those from the CT. The nodal coordinates obtained from the finite element model were used to find the identical locations on the referent CT to calculate the corresponding bone mineral density (BMD) values.

Results and Discussion

Piezoelectric Bone Remodeling Simulation

Starting from uniform density distribution, Figures 5A–C shows the evolution of the bone density at the intermediate and final time instants calculated using an interpolation post-processing technique. Regarding qualitative analysis, the final density distribution (see Figure 5C) was visually compared to the CT scan of the human proximal tibia (Figure 5D). The end configuration predicted fairly accurate density distribution with the trabecular bone beneath the tibial plateau, greater bone density in the medial region than in the lateral region, intramedullary canal with little trabecular bone, and cortical layers in the distal tibia. These observations are in good accordance with previous studies on tibial bone remodeling performed using commercial software (Rakotomanana, 2000; Perez et al., 2010; Robalo, 2011; Fang et al., 2013; Quilez et al., 2017). This provides a preliminary validation of the simulations performed with the open-source software framework illustrated in Figure 2A.

FIGURE 5

FIGURE 5. Evolution of bone density at time: (A)t = 100 days, (B)t = 200 days, and (C)t = 300 days. (D) Qualitative comparison of predicted bone density with the CT scan of the proximal tibia. (E) Location of the manually selected 182 nodal points for computing the RMS error and MD. (F) Quantitative comparison of predicted bone density with the BMD calculated from the CT (G) Histogram of absolute error (g/cm³) for the selected nodal points.

For quantitative analysis, RMS error and MD were computed for 182 nodal points (Figure 5E) that were selected using a non-probability convenience sampling technique (Lombardo et al., 2017). The HU values (−548 to 1883 HU) resulting from CT scans are in good accordance with those reported by Perez et al. (2007), Perez et al. (2010). Cortical bone was mostly observed in the cortex region, while trabecular bone was seen in the proximal tibia and only little in the medullary canal. The BMD calculated using the HU-density relationship described in Conversion of Hounsfield units into bone mineral density is shown in Figure 4F. By comparing the predicted (Figure 5C) and CT bone densities (Figure 5F) only for tibial domain (i.e., no side-plate nodes were considered), the RMS and MD values were computed to be 0.173 g/cm³ and 0.081 g/cm³, respectively. These quantitative comparisons provide supplementary validation of the implemented bone remodeling algorithm in open-source software. Further, for the selected nodal points, the histogram of the absolute errors is plotted in Figure 5G. In order to achieve better correspondence between simulation results and CT data, patient-specific distribution of uniform initial bone density, three-dimensional model with realistic boundary conditions, or a mixture of these should be considered. The findings presented here are intended to be representative of the tibial remodeling in general. Therefore, the approach used in this study is competent to predict the bone density distribution in 2D proximal tibia reasonably well.

Influence of the Uniform Initial Bone Density

Many bone remodeling models using different values of uniform initial bone density yield comparable but not identical final bone density distribution. In order to investigate the effect of varying initial bone density ${\rho }_0$ on the final density distribution, the model was analyzed for five different values of 0.2 g/cm³, 0.5 g/cm³, 0.8 g/cm³, 1.1 g/cm³, 1.4 g/cm³ and the obtained results are shown in Figure 6. As seen here, the final density distributions are not entirely identical and extensively rely on the value of uniform initial density. These findings are analogous to previous studies (Weinans et al., 1989; Orr, 1992; Turner et al., 1993) that also used the element-based remodeling algorithm.

FIGURE 6

FIGURE 6. Final bone density distributions (t = 300 days) obtained for distinct values of uniform initial density: (A)ρ₀= 0.2 g/cm³. (B)ρ₀= 0.5 g/cm³, (C)ρ₀= 0.8 g/cm³ (local regions of interests (ROIs) highlighted with white ellipses), (D)ρ₀= 1.1 g/cm³, and (E)ρ₀= 1.4 g/cm³.

At time 0 days, for the initial bone density equal to or higher than 0.8 g/cm³, the average bone density gradually declined throughout the remodeling period (see Figure 7A). On the contrary, when the same was equal to or lower than 0.5 g/cm³, the average bone density increased steadily in the first 75 days and then flattened out. Starting from different initial bone density, the resulting average bone densities tend to converge as time advances. On the 0th, 150th, and 300th day, the average bone density ranged from 0.2 to 1.4 g/cm³, 0.51–0.98 g/cm³, and 0.47–0.90 g/cm³, respectively. Further, the relative difference between these density values on the 300th day lies between 5.7 and 20.9%. Moreover, the bone density distributions obtained for different initial density values are also dissimilar (see Figure 6). In this figure, the region with the highest bone density (red) was reduced with an increase in initial density resulting in lower average bone density and this is in sync with the tendency noticed in Figure 7A.

FIGURE 7

FIGURE 7. Evolution of average bone density for: (A) distinct values of initial density ranging from 0.2 to 1.4 g/cm³ and (B) small variations within a range of ±1%, ±3%, and ±5% introduced to the selected initial bone density of 0.8 g/cm³.

As seen in Figure 7B, when variations within a range of ±1%, ±3%, and ±5% were introduced into the selected initial bone density of 0.8 g/cm³, the obtained density distribution at the final time was similar but not identical. It was thus observed that the bone density distribution (Figure 6) and the average bone density (Figures 7A,B) are dependent on the chosen value of initial density. However, the final density distribution shouldn’t be dependent on the initial conditions because they are mere numerical assumptions and not supported by substantial evidence. Therefore, to reduce such dependency, the lazy or dead zone (Carter, 1984; Huiskes et al., 1987) together with the saturated density change rate (Martínez-Reina et al., 2016; Su et al., 2019) could be incorporated into the remodeling algorithm implemented here.

In order to investigate the local bone adaptation, the evolution of average bone density in six regions of interest (ROIs) highlighted with white ellipses (see Figure 6C) was plotted in Figure 8. With the reference to the initial bone density of 0.8 g/cm³ in ROIs I, IV, and V, the average bone density increased for the lower values, whereas for the higher ones, it decreased over the remodeling period (see Figures 8A,D,E). Regarding ROIs I-IV, as days proceed, the differences in the average bone density were getting smaller until they were negligible at the final time (see Figures 8A–D). Also, a tendency to converge to a common value is remarkable. In the case of ROIs V and VI, the resulting average bone density deviated from each other over the bone remodeling span (see Figure 8 E and F). Therefore, this study shows that the final bone density distribution is location-specific. Here, the simulation results are based on the uniform initial bone density assumption opposing to the reality, but it would be meaningful for the study of numerical algorithms. Therefore, for interpreting the results, it is more logical to focus on the final bone density distribution instead of the evolution of density over time.

FIGURE 8

FIGURE 8. Average bone density evolution in the local regions of interest marked in Figure 6C: (A) ROI I (B) ROI II (C) ROI III (D) ROI IV (E) ROI V and (F) ROI VI.

Influence of the Reference Stimulus

The reference stimulus S_r has a considerable effect on bone remodeling simulations. The results obtained for S_r = 0.0001 J/g, 0.002 J/g, 0.004 J/g, 0.006 J/g, 0.008 J/g, and 0.01 J/g are depicted in Figure 9 A–F. It was observed that with an increase in the reference stimulus value, the highest bone density region (red) was reduced, and concomitantly the lowest bone density region (blue) was enlarged. Similar observations have been reported by Sarikanat and Yildiz (2011). In the present study, the reference stimulus of 0.004 J/g was selected as it provides the density distribution, which is in good agreement with the real bone. For the reference stimulus values less than 0.004 J/g, the average bone densities increased steadily, while for greater values, the average bone densities decreased steadily at a similar rate (Figure 9G). On the 300th day, for the reference stimulus values from 0.0001 J/g to 0.01 J/g, the resulting average bone density ranges from 1.12 g/cm³ to 0.28 g/cm³. It is remarkable that a lower value of the reference stimulus resulted in higher average bone density. Therefore, the bone density distribution and the average bone density (Figure 9) are greatly dependent on the reference stimulus value. However, in reality, bone tissue is not very sensitive to changes in the reference mechanical stimulus. Taking this into account, many studies have considered a nonlinear remodeling rate relation including the lazy or dead zone (Carter, 1984; Huiskes et al., 1987; Fischer et al., 1995).

FIGURE 9

FIGURE 9. Final bone density distributions (t = 300 days) obtained for distinct values of reference stimulus: (A)S_r = 0.0001 J/g, (B)S_r = 0.002 J/g (C)S_r = 0.004 J/g (selected), (D)S_r = 0.006 J/g (E)S_r = 0.008 J/g, and (F)S_r = 0.010 J/g. (G) Evolution of average bone density obtained for the reference stimulus values ranging from 0.0001 J/g to 0.01 J/g.

Bone Piezoelectricity

It is evident from experimental studies (Fukada and Yasuda, 1957; Bassett and Becker, 1962; Becker, 1964; Bassett et al., 1964; Black and Korostoff, 1974; Yasuda, 1977; Nade, 1994; Zigman et al., 2013) that the mechanical loading induces changes in the bone electric potentials in a way that regions exposed to compressive loads generated negative potentials, whereas those exposed to tensile loads generated positive potentials. Figure 10 demonstrates changes in the distribution and magnitude of the electric potentials generated at different time instants of the day due to the application of varying mechanical loads, representing the direct piezoelectric effect. For better comparison within the generated electrical potentials, they were normalized with respect to their mean amplitude. For electromechanical simulations, negative potentials are associated with osteoblast-induced bone formation, whereas positive potentials are associated with osteoclast-induced bone resorption (Qin and Ye, 2004; Fernández et al., 2012a; Cerrolaza et al., 2017). These electrical potentials play a vital role in the process of bone healing and remodeling (Fukada and Yasuda, 1957; Fukada and Yasuda, 1964; Marino and Becker, 1970; Becker and Spadaro, 1972; Gjelsvik, 1973; Ramtani, 2008; Zigman et al., 2013). A sensitivity analysis was conducted to investigate the influence of piezoelectric and permittivity tensors on the generated electric potentials. This analysis demonstrated that the model was very sensitive to the values of these parameters and similar findings have been observed by Fernández et al. (2012c).

FIGURE 10

FIGURE 10. Normalized electric potentials generated at various time instants of the 238th day due to the application of varying mechanical loads during walking.

Therapeutic Electrical Stimulation to Reduce Bone Loss

The bone density distribution predicted after the remodeling period (Figure 5C) was assumed to be an initial state for the simulations of reduced physical activity and therapeutic electrical stimulation (see Loads and boundary conditions for more details). The density distribution predicted after the reduced physical activity period is shown in Figure 11A. In order to investigate the local variations in bone density distributions as a result of the reduced physical activity, the differences between bone density after the remodeling period (i.e., on 300th day) and after the reduced physical activity period (i.e., on 400th day) are plotted in Figure 11B. It is notable that the density distribution in both condyles and the intramedullary canal was lower than that obtained after the remodeling period.

FIGURE 11

FIGURE 11. (A) Predicted bone density distribution after the reduced physical activity period (t = 400 days). (B) Differences in predicted bone density, after the remodeling period (t = 300 days) and after the reduced physical activity period (t = 400 days). (C) Electric potential generated at time t = 400 days, when a surface electric charge was applied to the medial condyle for the reduced physical activity period (i.e., 300–400 days). (D) Differences in predicted bone density at time t = 400 days, when therapeutic electrical stimulation was applied along with acting mechanical loads and when only mechanical loads were acting.

The most significant benefit of including the piezoelectricity in bone remodeling study is the ability to change bone density under electrical stimulation (Fernández et al., 2012a) through the converse piezoelectric effect. When a surface electric charge was applied to the tibia during the reduced physical activity period, a negative potential of nearly −36 V was noticed (see Figure 11C). In order to show the resulting change in bone density, differences between densities while the electric charge was applied in addition to the acting mechanical loads and while only the mechanical loads were acting, are plotted in Figure 11D. Minor variations in bone density were observed, where areas in red and blue correspond to growth and reduction in bone density, respectively. In comparison with the reduction in bone density in condyles and intramedullary canal owing to reduced physical activity (Figure 11B), the therapeutic electrical stimulation affects only the area in its vicinity, resulting in a clear increase in the density between the condyles (Figure 11D). These preliminary observations showed that the electromechanical loads affect the bone density evolution. Hence, therapeutic electrical stimulation could be treated as a supplement to improve bone remodeling and to minimize bone loss, e.g., physical impairment or osteoporosis. Direct comparison of the simulation predictions with the experimental data is not feasible because of its unavailability. Matrix piezoelectricity and streaming potential are the two mechanisms responsible for the piezoelectricity in bone (Fernández et al., 2012a; Cerrolaza et al., 2017; Mohammedkhah et al., 2019). However, in the present study, small changes in bone density after electrical stimulation were ascribed only to the matrix piezoelectricity. Therefore, the obtained results inspire further advancement of a multi-physics model that involves multiple coupled phenomena such as matrix piezoelectricity (Fukada and Yasuda, 1957), strain-generated fluid flow (Cowin et al., 1995), and streaming potential (Pollack et al., 1984) leads to more promising results for therapeutic electrical stimulation.

So far, only a small number of bone remodeling models have been studied considering the piezoelectric effect and that even using either commercial software or in-house programs that are not easily accessible. Nonetheless, in the current study, the piezoelectric strain-adaptive bone remodeling in the human proximal tibia was analyzed for the first time and that too using an open-source framework. The Python code files developed are available on GitHub for easy access. Here, the most contemporary model of bone remodeling that considers the piezoelectric properties of bone was implemented to investigate bone remodeling under electrical stimulation. The predicted bone density distributions were validated qualitatively by visually comparing with the radiographic scan and quantitatively by calculating the RMS error between the predicted BMD and the BMD obtained from the CT. Both direct and converse piezoelectric effects were shown by combining the effects of mechanical and electric fields on the response of bone. A detailed parametric study was carried out to investigate the influence of reference stimulus and uniform initial bone density on the final density distribution. For time discretization in bone remodeling studies, the explicit Euler method was commonly used, which is not very precise and for large time-step size, becomes unstable. Thus, in this study, a more stable and accurate fourth-order Runge-Kutta method was used. The findings of this study highlight the role of bone piezoelectricity in the therapeutic electrical stimulation and could be used clinically to improve the osseointegration of the electrically active implants.

Nonetheless, this study has some drawbacks. In spite of the substantial contribution of matrix piezoelectricity and streaming potential to the electromechanical properties of bone (Fernández et al., 2012c), their coupling is not yet well understood. In the literature, there are many discrepancies in the measured values of the piezoelectric strain tensors (Fukada, 1968; Fotiadis et al., 1999; Qin and Ye, 2004; Mohammadkhah et al., 2019). Here, for computational modeling purposes, the bone was assumed to be electrically homogeneous, which is not the case in reality (Johnson et al., 1980). This study does not take into account the external bone remodeling (Goda et al., 2016). The implemented algorithm does not consider the effect of gender, age, disease, and injury on the bone remodeling process. The findings of the present study cannot be generalized to the severely osteoporotic or fractured bone (i.e., with very less or no physical activity), as the piezoelectric bone remodeling model suggests replacing only the part of mechanical load resulting from the daily physical activity by electrical stimulation.

In future work, the constitutive laws for mechanical and electrical behavior of bone will be adapted to include coupled multi-physics phenomena such as matrix piezoelectricity and streaming potential to better predict bone density evolution. The simulations of reduced physical activity and electrical stimulation can be performed considering the initial bone density distribution to be patient-specific. The electric stimulation of varying magnitude will be applied to the different regions of the tibia investigating the effect of stimulation intensity and site of excitation on bone remodeling. To correlate HU values with BMD, a calibration phantom could be used instead of a calibration function. By including the results from musculoskeletal multibody simulations, the boundary conditions should be improved to perform full 3D simulations. From a rehabilitation perspective, future studies should investigate whether the numerical model is indicative of the potential therapeutic values of electrical stimulation. The implemented piezoelectric bone remodeling algorithm can also be employed for applications investigating the effect of electrically active implants (Schmidt et al., 2015; Zimmermann and van Rienen, 2015; Raben et al., 2019; Zimmermann et al., 2021) in the adjacent bone tissue with respect to peri-implant bone remodeling.

Conclusion

This is the first attempt to investigate the piezoelectric bone remodeling in the human tibia for daily walking activity using an open-source framework. The simulation results predicted reasonably accurate bone density distributions that were validated both qualitatively and quantitatively against the CT data. The parametric analysis showed that different uniform initial bone density and reference stimulus values resulted in different density distribution at the final time. A reduction in bone density was observed for the reduced physical activity compared to the daily physical activity. Therapeutic electrical stimulation applied over a period of reduced physical activity showed increased bone deposition suggesting that in the case of bone loss, e.g., physical impairment or osteoporosis, the mechanical loads can be replaced in part by electrical stimuli that enhance bone density or reduce bone loss.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/YDBansod/Bone_Remodelling.

Ethics Statement

The studies involving human participants were reviewed and approved by Ethics Committee at the Medical Faculty of the University of Rostock, Germany. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

Author Contributions

UvR and RB were responsible for resources and funding acquisition. YDB and MK have prepared the initial draft of the manuscript. YDB and MK have performed the simulations. All authors contributed to the analysis of the results. DK, RB, and UvR have critically revised the manuscript for important suggestions. All authors have read and approved the final manuscript to be submitted.

Funding

This study was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) SFB 1270/1 ELAINE-299150580.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors gratefully acknowledge support from the German Research Foundation (DFG) via the Collaborative Research Center 1270 ELAINE—Electrically Active Implants. The DFG had no role in the study design, data collection, data analysis, decision to publish, and preparation of the manuscript.

References

Ahrens, J., Geveci, B., and Law, C. (2005). ParaView: An End-User Tool for Large-Data Visualization. The Visualization Handbook 717, 717–731. doi:10.1016/B978-012387582-2/50038-1

CrossRef Full Text | Google Scholar

Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., et al. (2015). The FEniCS Project Version 1.5. Archive Numer. Softw. 3, 9–23. doi:10.11588/ans.2015.100.20553

CrossRef Full Text | Google Scholar

Bansod, Y. D., Kebbach, M., Kluess, D., Bader, R., and van Rienen, U. (2021). Finite Element Analysis of Bone Remodelling with Piezoelectric Effects Using an Open-Source Framework. Biomech. Model. Mechanobiol. 20, 1147–1166. doi:10.1007/s10237-021-01439-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Bansod, Y. D., and van Rienen, U. (2019). “Numerical Analysis of Electromechanically Driven Bone Remodeling Using the Open-Source Software Framework,” in 41st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), 23-27 July 2019 (Berlin, Germany: IEEE), 6466–6471. doi:10.1109/EMBC.2019.8856543

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassett, C. A. L., and Becker, R. O. (1962). Generation of Electric Potentials by Bone in Response to Mechanical Stress. Science 137, 1063–1064. doi:10.1126/science.137.3535.1063

PubMed Abstract | CrossRef Full Text | Google Scholar

Bassett, C. A. L., Pawluk, R. J., and Becker, R. O. (1964). Effects of Electric Currents on Bone In Vivo. Nature 204, 652–654. doi:10.1038/204652a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Batra, R. C., and Yang, J. S. (1995). Saint-Venant's Principle in Linear Piezoelectricity. J. Elasticity 38, 209–218. doi:10.1007/BF00042498

CrossRef Full Text | Google Scholar

Becker, R. O., and Spadaro, J. A. (1972). Electrical Stimulation of Partial Limb Regeneration in Mammals. Bull. N. Y Acad. Med. 48, 627–641.

PubMed Abstract | Google Scholar

Becker, R. O. (1964). Bioelectric Factors Controlling Bone Structure. Bone Biodynamics.

Black, J., and Korostoff, E. (1974). Strain-Related Potentials in Living Bone. Ann. NY Acad. Sci. 238, 95–120. doi:10.1111/j.1749-6632.1974.tb26781.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bouguecha, A., Weigel, N., Behrens, B.-A., Stukenborg-Colsman, C., and Waizy, H. (2011). Numerical Simulation of Strain-Adaptive Bone Remodelling in the Ankle Joint. BioMedical Eng. OnLine 10 (1), 58–13. doi:10.1186/1475-925X-10-58

CrossRef Full Text | Google Scholar

Carter, D. R. (1984). Mechanical Loading Histories and Cortical Bone Remodeling. Calcif Tissue Int. 36, S19–S24. doi:10.1007/BF02406129

PubMed Abstract | CrossRef Full Text | Google Scholar

Cerrolaza, M., Duarte, V., and Garzón-Alvarado, D. (2017). Analysis of Bone Remodeling under Piezoelectricity Effects Using Boundary Elements. J. Bionic. Eng. 14, 659–671. doi:10.1016/S1672-6529(16)60432-8

CrossRef Full Text | Google Scholar

Cowin, S. (2001). Bone Mechanics Handbook. Boca Raton: CRC Press. doi:10.1201/b14263

CrossRef Full Text

Cowin, S. C., Weinbaum, S., and Zeng, Y. (1995). A Case for Bone Canaliculi as the Anatomical Site of Strain Generated Potentials. J. Biomech. 28, 1281–1297. doi:10.1016/0021-9290(95)00058-P

CrossRef Full Text | Google Scholar

Duarte, V. d. J., Thoeni, K., Garzón-Alvarado, D., and Cerrolaza, M. (2018). A Simplified Scheme for Piezoelectric Anisotropic Analysis in Human Vertebrae Using Integral Methods. Math. Probl. Eng. 2018, 1–8. doi:10.1155/2018/2783867

CrossRef Full Text | Google Scholar

Duda, G. N., Mandruzzato, F., Heller, M., Goldhahn, J., Moser, R., Hehli, M., et al. (2001). Mechanical Boundary Conditions of Fracture Healing: Borderline Indications in the Treatment of Unreamed Tibial Nailing. J. Biomech. 34, 639–650. doi:10.1016/S0021-9290(00)00237-2

CrossRef Full Text | Google Scholar

Fang, J., Gong Kong Zhu, H. L. D., Kong, L., and Zhu, D. (2013). Simulation on the Internal Structure of Three-Dimensional Proximal Tibia under Different Mechanical Environments. BioMedical Eng. OnLine 12, 1–17. doi:10.1186/1475-925X-12-130

CrossRef Full Text | Google Scholar

Fernández, J. R., García-aznar, J. M., and Martínez, R. (2012a). Numerical Analysis of a Piezoelectric Bone Remodelling Problem. Eur. J. Appl. Math. 23, 635–657. doi:10.1017/S0956792512000150

CrossRef Full Text | Google Scholar

Fernández, J. R., García-Aznar, J. M., and Martínez, R. (2012c). Piezoelectricity Could Predict Sites of Formation/Resorption in Bone Remodelling and Modelling. J. Theor. Biol. 292, 86–92. doi:10.1016/j.jtbi.2011.09.032

CrossRef Full Text | Google Scholar

Fernández, J. R., García-Aznar, J. M., Martínez, R., and Viaño, J. M. (2010). Numerical Analysis of a Strain-Adaptive Bone Remodelling Problem. Comp. Methods Appl. Mech. Eng. 199, 1549–1557. doi:10.1016/j.cma.2010.01.005

CrossRef Full Text | Google Scholar

Fernández, R. M. (2010). “Numerical Analysis and Simulations in Bone Remodeling Models,”. PhD dissertation (Spain: Universidade de Santiago de Compostela).

Google Scholar

Fischer, K. J., Jacobs, C. R., and Carter, D. R. (1995). Computational Method for Determination of Bone and Joint Loads Using Bone Density Distributions. J. Biomech. 28, 1127–1135. doi:10.1016/0021-9290(94)00182-4

CrossRef Full Text | Google Scholar

Fischer, K. J., Jacobs, C. R., Levenston, M. E., and Carter, D. R. (1997). Observations of Convergence and Uniqueness of Node-Based Bone Remodeling Simulations. Ann. Biomed. Eng. 25, 261–268. doi:10.1007/BF02648040

PubMed Abstract | CrossRef Full Text | Google Scholar

Fotiadis, D. I., Foutsitzi, G., and Massalas, C. V. (1999). Wave Propagation Modeling in Human Long Bones. Acta Mechanica 137, 65–81. doi:10.1007/BF01313145

CrossRef Full Text | Google Scholar

Fukada, E. (1968). Mechanical Deformation and Electrical Polarization in Biological Substances. Bir. 5, 199–208. doi:10.3233/BIR-1968-5302

CrossRef Full Text | Google Scholar

Fukada, E., and Yasuda, I. (1957). On the Piezoelectric Effect of Bone. J. Phys. Soc. Jpn. 12, 1158–1162. doi:10.1143/JPSJ.12.1158

CrossRef Full Text | Google Scholar

Fukada, E., and Yasuda, I. (1964). Piezoelectric Effects in Collagen. Jpn. J. Appl. Phys. 3, 117–121. doi:10.1143/jjap.3.117

CrossRef Full Text | Google Scholar

García-Aznar, J. M., Rueberg, T., and Doblare, M. (2005). A Bone Remodelling Model Coupling Microdamage Growth and Repair by 3D BMU-Activity. Biomech. Model. Mechanobiol. 4, 147–167. doi:10.1007/s10237-005-0067-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Garijo, N., Fernández, J. R., Pérez, M. A., and García-Aznar, J. M. (2014). Numerical Stability and Convergence Analysis of Bone Remodeling Model. Comp. Methods Appl. Mech. Eng. 271, 253–268. doi:10.1016/j.cma.2013.12.014

CrossRef Full Text | Google Scholar

Garijo, N., Verdonschot, N., Engelborghs, K., García-Aznar, J. M., and Pérez, M. A. (2017). Subject-Specific Musculoskeletal Loading of the Tibia: Computational Load Estimation. J. Mech. Behav. Biomed. Mater. 65, 334–343. doi:10.1016/j.jmbbm.2016.08.026

CrossRef Full Text | Google Scholar

Garzón-Alvarado, D. A., Ramírez-Martínez, A. M., and Cardozo de Martínez, C. (2012). Numerical Test Concerning Bone Mass Apposition under Electrical and Mechanical Stimulus. Theor. Biol. Med. Model. 9, 14. doi:10.1186/1742-4682-9-14

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerhard, F. A., Webster, D. J., Van Lenthe, G. H., and Müller, R. (2009). In Silico Biology of Bone Modelling and Remodelling: Adaptation. Phil. Trans. R. Soc. A. 367, 2011–2030. doi:10.1098/rsta.2008.0297

PubMed Abstract | CrossRef Full Text | Google Scholar

Geuzaine, C., and Remacle, J.-F. (2009). Gmsh: A 3-D Finite Element Mesh Generator with Built-In Pre- and post-processing Facilities. Int. J. Numer. Meth. Engng. 79, 1309–1331. doi:10.1002/nme.2579

CrossRef Full Text | Google Scholar

Gjelsvik, A. (1973). Bone Remodeling and Piezoelectricity - I. J. Biomech. 6, 69–77. doi:10.1016/0021-9290(73)90039-0

CrossRef Full Text | Google Scholar

Goda, I., Ganghoffer, J.-F., and Maurice, G. (2016). Combined Bone Internal and External Remodeling Based on Eshelby Stress. Int. J. Sol. Structures 94-95, 138–157. doi:10.1016/j.ijsolstr.2016.04.036

CrossRef Full Text | Google Scholar

Gong, H., Kong, L., Zhang, R., Fang, J., and Zhao, M. (2013). A Femur-Implant Model for the Prediction of Bone Remodeling Behavior Induced by Cementless Stem. J. Bionic. Eng. 10, 350–358. doi:10.1016/S1672-6529(13)60230-9

CrossRef Full Text | Google Scholar

González-Carbonell, R. A., Ortiz-Prado, A., Cisneros-Hidalgo, Y. A., and Alpizar-Aguirre, A. (2015). “Bone Remodeling Simulation of Subject-specific Model of Tibia under Torque,” in VI Latin American Congress on Biomedical Engineering CLAIB 2014, Paraná, Argentina. October 2014 (Cham: Springer), 305–308. doi:10.1007/978-3-319-13117-7_79

CrossRef Full Text | Google Scholar

Hambli, R., Katerchi, H., and Benhamou, C.-L. (2011). Multiscale Methodology for Bone Remodelling Simulation Using Coupled Finite Element and Neural Network Computation. Biomech. Model. Mechanobiol. 10, 133–145. doi:10.1007/s10237-010-0222-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hastings, G. W., and Mahmud, F. A. (1988). Electrical Effects in Bone. J. Biomed. Eng. 10, 515–521. doi:10.1016/0141-5425(88)90109-4

CrossRef Full Text | Google Scholar

Huiskes, R., Weinans, H., Grootenboer, H. J., Dalstra, M., Fudala, B., and Slooff, T. J. (1987). Adaptive Bone-Remodeling Theory Applied to Prosthetic-Design Analysis. J. Biomech. 20, 1135–1150. doi:10.1016/0021-9290(87)90030-3

CrossRef Full Text | Google Scholar

Johnson, M. W., Williams, W. S., and Gross, D. (1980). Ceramic Models for Piezoelectricity in Dry Bone. J. Biomech. 13, 565–573. doi:10.1016/0021-9290(80)90057-3

CrossRef Full Text | Google Scholar

Knowles, N. K., Reeves, J. M., and Ferreira, L. M. (2016). Quantitative Computed Tomography (Qct) Derived Bone Mineral Density (BMD) in Finite Element Studies: A Review of the Literature. J. Exp. Ortop. 3, 1–16. doi:10.1186/s40634-016-0072-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Lian, Q., Li, D., Jin, Z., Wang, Z., and Sun, Y. (2014). Patient-Specific Design and Biomechanical Evaluation of a Novel Bipolar Femoral Hemi-Knee Prosthesis. J. Bionic. Eng. 11, 259–267. doi:10.1016/S1672-6529(14)60039-1

CrossRef Full Text | Google Scholar

Logg, A., Mardal, K. A., and Wells, G. (2012). Automated Solution of Differential Equations by the Finite Element Method: The Fenics Book. Berlin Heidelberg: Springer Science & Business Media. doi:10.1007/978-3-642-23099-8

CrossRef Full Text

Lombardo, G., Pighi, J., Marincola, M., Corrocher, G., Simancas-PallaresNocini, M. P. F., and Nocini, P. F. (2017). Cumulative Success Rate of Short and Ultrashort Implants Supporting Single Crowns in the Posterior Maxilla: A 3-Year Retrospective Study. Int. J. Dentistry 2017, 1–10. doi:10.1155/2017/8434281

PubMed Abstract | CrossRef Full Text | Google Scholar

Marino, A. A., and Becker, R. O. (1970). Piezoelectric Effect and Growth Control in Bone. Nature 228, 473–474. doi:10.1038/228473a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Martínez-Reina, J., Ojeda, J., and Mayo, J. (2016). On the Use of Bone Remodelling Models to Estimate the Density Distribution of Bones. Uniqueness of the Solution. Plos One 11, e0148603. doi:10.1371/journal.pone.0148603

PubMed Abstract | CrossRef Full Text | Google Scholar

Mauck, J., Wieding, J., Kluess, D., and Bader, R. (2016). Numerical Simulation of Mechanically Stimulated Bone Remodelling. Curr. Dir. Biomed. Eng. 2, 643–647. doi:10.1515/cdbme-2016-0141

CrossRef Full Text | Google Scholar

Minary-Jolandan, M., and Yu, M.-F. (2009). Nanoscale Characterization of Isolated Individual Type I Collagen Fibrils: Polarization and Piezoelectricity. Nanotechnology 20, 085706. doi:10.1088/0957-4484/20/8/085706

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohammadkhah, M., Marinkovic, D., Zehn, M., and Checa, S. (2019). A Review on Computer Modeling of Bone Piezoelectricity and its Application to Bone Adaptation and Regeneration. Bone 127, 544–555. doi:10.1016/j.bone.2019.07.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Nade, S. (1994). Stimulating Osteogenesis. Injury 25, 577–583. doi:10.1016/0020-1383(94)90031-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Nyman, J. S., Hazelwood, S. J., Rodrigo, J. J., Martin, R. B., and Yeh, O. C. (2004). Long Stemmed Total Knee Arthroplasty with Interlocking Screws: A Computational Bone Adaptation Study. J. Orthop. Res. 22, 51–57. doi:10.1016/S0736-0266(03)00159-1

CrossRef Full Text | Google Scholar

Orr, T. E. (1992). “The Role of Mechanical Stresses in Bone Remodelling,”. PhD dissertation (Stanford: Stanford University).

Google Scholar

Papathanasopoulou, V. A., Fotiadis, D. I., Foutsitzi, G., and Massalas, C. V. (2002). A Poroelastic Bone Model for Internal Remodeling. Int. J. Eng. Sci. 40, 511–530. doi:10.1016/S0020-7225(01)00076-3

CrossRef Full Text | Google Scholar

Peng, L., Bai, J., Zeng, X., and Zhou, Y. (2006). Comparison of Isotropic and Orthotropic Material Property Assignments on Femoral Finite Element Models under Two Loading Conditions. Med. Eng. Phys. 28, 227–233. doi:10.1016/j.medengphy.2005.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Pérez, M. A., Fornells, P., Doblaré, M., and García-Aznar, J. M. (2010). Comparative Analysis of Bone Remodelling Models with Respect to Computerised Tomography-Based Finite Element Models of Bone. Comp. Methods Biomech. Biomed. Eng. 13, 71–80. doi:10.1080/10255840903045029

CrossRef Full Text | Google Scholar

Perez, M. A., Fornells, P., Garcıa-Aznar, J. M., and Doblare, M. (2007). Validation of Bone Remodelling Models Applied to Different Bone Types Using Mimics. Zaragoza, Spain: University of Zaragoza.

Pienkowski, D., and Pollack, S. R. (1983). The Origin of Stress-Generated Potentials in Fluid-Saturated Bone. J. Orthop. Res. 1, 30–41. doi:10.1002/jor.1100010105

CrossRef Full Text | Google Scholar

Pollack, S. R., Petrov, N., Salzstein, R., Brankov, G., and Blagoeva, R. (1984). An Anatomical Model for Streaming Potentials in Osteons. J. Biomech. 17, 627–636. doi:10.1016/0021-9290(84)90094-0

CrossRef Full Text | Google Scholar

Qin, Q.-H. (2003). Fracture Analysis of Cracked Thermopiezoelectric Materials by BEM. Ejbe 1, 283–301. doi:10.14713/ejbe.v1i2.762

CrossRef Full Text | Google Scholar

Qin, Q.-H., Qu, C., and Ye, J. (2005). Thermoelectroelastic Solutions for Surface Bone Remodeling under Axial and Transverse Loads. Biomaterials 26, 6798–6810. doi:10.1016/j.biomaterials.2005.03.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, Q.-H., and Ye, J.-Q. (2004). Thermoelectroelastic Solutions for Internal Bone Remodeling under Axial and Transverse Loads. Int. J. Sol. Structures 41, 2447–2460. doi:10.1016/j.ijsolstr.2003.12.026

CrossRef Full Text | Google Scholar

Qu, C.-Y., and Yu, S.-W. (2011). The Damage and Healing of Bone in the Disuse State under Mechanical and Electro-Magnetic Loadings. Proced. Eng. 10, 171–176. doi:10.1016/j.proeng.2011.04.031

CrossRef Full Text | Google Scholar

Qu, C., Qin, Q.-H., and Kang, Y. (2006). A Hypothetical Mechanism of Bone Remodeling and Modeling under Electromagnetic Loads. Biomaterials 27, 4050–4057. doi:10.1016/j.biomaterials.2006.03.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Quilez, M. P., Seral, B., and Pérez, M. A. (2017). Biomechanical Evaluation of Tibial Bone Adaptation after Revision Total Knee Arthroplasty: A Comparison of Different Implant Systems. Plos one 12, e0184361. doi:10.1371/journal.pone.0184361

PubMed Abstract | CrossRef Full Text | Google Scholar

Raben, H., Kämmerer, P. W., Bader, R., and van Rienen, U. (2019). Establishment of a Numerical Model to Design an Electro-Stimulating System for a Porcine Mandibular Critical Size Defect. Appl. Sci. 9, 2160. doi:10.3390/app9102160

CrossRef Full Text | Google Scholar

Rakotomanana, L. (2000). “The Long-Term Behaviour of Human Joints with Orthopaedic Prosthesis: Finite Element Models,” in European Congress on Comput. Meth. In Appl. Sci. Eng. (Barcelona, Spain: ECCOMAS).

Google Scholar

Ramtani, S. (2008). Electro-Mechanics of Bone Remodeling. Int. J. Eng. Sci. 46, 1173–1182. doi:10.1016/j.ijengsci.2008.06.001

CrossRef Full Text | Google Scholar

Ribes, A., and Caremoli, C. (2007). “Salome Platform Component Model for Numerical Simulation,” in 31st Annual International Computer Software and Applications Conference (COMPSAC 2007), 24-27 July 2007 (Beijing: IEEE), 553–564. doi:10.1109/COMPSAC.2007.185qx

CrossRef Full Text | Google Scholar

Robalo, T. (2011). Analysis of Bone Remodeling in the Tibia after Total Knee Prosthesis. Lisbon, Portugal: Instituto Superior Técnico, Universida de Técnica de Lisboa.

Robling, A. G., and Turner, C. H. (2009). Mechanical Signaling for Bone Modeling and Remodeling. Crit. Rev. Eukar Gene Expr. 19, 319–338. doi:10.1615/CritRevEukarGeneExpr.v19.i4.50

CrossRef Full Text | Google Scholar

Sarikanat, M., and Yildiz, H. (2011). Determination of Bone Density Distribution in Proximal Femur by Using the 3D Orthotropic Bone Adaptation Model. Proc. Inst. Mech. Eng. H 225, 365–375. doi:10.1177/09544119JEIM895

CrossRef Full Text | Google Scholar

Schmidt, C., Zimmermann, U., and van Rienen, U. (2015). Modeling of an Optimized Electrostimulative Hip Revision System under Consideration of Uncertainty in the Conductivity of Bone Tissue. IEEE J. Biomed. Health Inform. 19, 1321–1330. doi:10.1109/JBHI.2015.2423705

CrossRef Full Text | Google Scholar

Shirazi Beheshtiha, A., and Nackenhorst, U. (2015). Computational Simulation of Piezo-Electrically Stimulated Bone Remodeling Surrounding Teeth Implant. Proc. Appl. Math. Mech. 15, 111–112. doi:10.1002/pamm.201510046

CrossRef Full Text | Google Scholar

Su, K., Yuan, L., Yang, J., and Du, J. (2019). Numerical Simulation of Mandible Bone Remodeling under Tooth Loading: A Parametric Study. Sci. Rep. 9, 1–12. doi:10.1038/s41598-019-51429-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, Y., Kluess, D., Mittelmeier, W., van Rienen, U., and Bader, R. (2016). An Automatic Approach for Calibrating Dielectric Bone Properties by Combining Finite-Element and Optimization Software Tools. Comp. Methods Biomech. Biomed. Eng. 19, 1306–1313. doi:10.1080/10255842.2015.1131980

CrossRef Full Text | Google Scholar

Turner, C. H., Anne, V., and Pidaparti, R. M. (1993). Convergence and Uniqueness of Bone Adaptation Models. Trans. Orthop. Res. Soc. 18, 527.

Google Scholar

Wang, X., Thomas, C. D. L., Clement, J. G., Das, R., Davies, H., and Fernandez, J. W. (2016). A Mechanostatistical Approach to Cortical Bone Remodelling: An Equine Model. Biomech. Model. Mechanobiol. 15, 29–42. doi:10.1007/s10237-015-0669-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Weinans, H., Huiskes, R., and Grootenboer, H. J. (1989). Convergence and Uniqueness of Adaptive Bone Remodeling. Trans. Orthop. Res. Soc. 310, 10003529282.

Google Scholar

Weinans, H., Huiskes, R., and Grootenboer, H. J. (1992). The Behavior of Adaptive Bone-Remodeling Simulation Models. J. Biomech. 25, 1425–1441. doi:10.1016/0021-9290(92)90056-7

CrossRef Full Text | Google Scholar

Wirtz, D. C., Schiffers, N., Pandorf, T., Radermacher, K., Weichert, D., and Forst, R. (2000). Critical Evaluation of Known Bone Material Properties to Realize Anisotropic FE-Simulation of the Proximal Femur. J. Biomech. 33, 1325–1330. doi:10.1016/S0021-9290(00)00069-5

CrossRef Full Text | Google Scholar

Wolff, J. (1893). Das Gesetz der Transformation der Knochen. Dtsch Med. Wochenschr 19, 1222–1224. doi:10.1055/s-0028-1144106

CrossRef Full Text | Google Scholar

Yasuda, I. (1954). On the Piezoelectric Activity of Bone. J. Jpn. Orthop. Surg. Soc. 28, 267.

Google Scholar

Yasuda, I. W. (1977). Electrical Callus and Callus Formation by Electret. Clin. Orthopaedics Relat. Res. 124, 53–56. doi:10.1097/00003086-197705000-00007

CrossRef Full Text | Google Scholar

Yushkevich, P. A., Piven, J., Hazlett, H. C., Smith, R. G., Ho, S., Gee, J. C., et al. (2006). User-Guided 3D Active Contour Segmentation of Anatomical Structures: Significantly Improved Efficiency and Reliability. Neuroimage 31, 1116–1128. doi:10.1016/j.neuroimage.2006.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Zigman, T., Davila, S., Dobric, I., Antoljak, T., Augustin, G., Rajacic, D., et al. (2013). Intraoperative Measurement of Bone Electrical Potential: A Piece in the Puzzle of Understanding Fracture Healing. Injury 44, S16–S19. doi:10.1016/S0020-1383(13)70191-8

CrossRef Full Text | Google Scholar

Zimmermann, U., Ebner, C., Su, Y., Bender, T., Bansod, Y. D., Mittelmeier, W., et al. (2021). Numerical Simulation of Electric Field Distribution Around an Instrumented Total Hip Stem. Appl. Sci. 11, 6677. doi:10.3390/app11156677

CrossRef Full Text | Google Scholar

Zimmermann, U., and van Rienen, U. (2015). “Investigation of the Electric Field Distribution within an Electro-Stimulated Microstructure of Cancellous Bone,” in International Conference on the Computation of Electromagnetic Fields (Compumag 2015), Montreal, Canada, 28 June - 2 July 2015, 475–476.

Google Scholar