Numerical test concerning bone mass apposition under electrical and mechanical stimulus

Abstract

This article proposes a model of bone remodeling that encompasses mechanical and electrical stimuli. The remodeling formulation proposed by Weinans and collaborators was used as the basis of this research, with a literature review allowing a constitutive model evaluating the permittivity of bone tissue to be developed. This allowed the mass distribution that depends on mechanical and electrical stimuli to be obtained. The remaining constants were established through numerical experimentation. The results demonstrate that mass distribution is altered under electrical stimulation, generally resulting in a greater deposition of mass. In addition, the frequency of application of an electric field can affect the distribution of mass; at a lower frequency there is more mass in the domain. These numerical experiments open up discussion concerning the importance of the electric field in the remodeling process and propose the quantification of their effects.

Background

Bones provide mechanical stability to the human body and are a source of minerals for metabolism [1]. Bones have been studied extensively from the mechanical and mineral standpoint, and in terms of functionality [1, 2]. From a mechanical point of view, bones can be adapted to loads on trajectories of stress through mineral apposition, which is due to the action of osteoblasts [1–4]. Furthermore, they reabsorb minerals when the mechanical stimulus is sufficiently low, as it is unnecessary to maintain structure [2]. Reabsorption is directed by osteoclasts. Osteoblasts and osteoclasts are the primary cells involved during bone remodeling that are stimulated by the action of mechanical strain sensors, for example, osteocytes [2]. These three cell types play an important role during the processes of replacement, maintenance and modeling of bones [1].

Following the work of Meyer during the nineteenth century, Wolff [5] proposed a theory of trabecular bone architecture, which assumes that trajectories of high mechanical stress form the trabecular bone. In 1987, Frost [6–8] suggested an adaptive mechanism of bone mass as a function of mechanical stress. Consequently, several bone remodeling algorithms have been developed including those proposed by Frost [8], Pauwels [9], Kummer [10], Cowin [11–13] and Hegedus [14], which predict the formation of bone structure from internal mechanical loads studied in terms of stress and strain.

From mechanical models of bone remodeling, sophisticated studies have been carried out concerning the processes of apposition and reabsorption during bone turnover, and particularly concerning the distribution of mass in the femur [15, 16], hip replacement [17, 18], and dental implants [19]. Generally, these studies have been phenomenological. Therefore, researchers have made significant efforts to include mathematical models and the role of cell biology and biochemistry in the remodeling process, resulting in research that begins at the microscopic level, concerning the effects of basic cellular remodeling units (BMU, Basic multicellular units) during tissue replacement [20, 21]. From the perspective of BMU, important work was initiated at the biochemical and mechanical level concerning the effects of cracks [22], cell cycles throughout adult life [23], active molecules within each cell [24] and the spatial distribution of each BMU [25]. With these important advances in the understanding of bone remodeling, researchers in the field increasingly turned to the study of other biophysical stimuli that can affect this process. Most models have not taken account of the physical-chemical phenomena of tissue mechano-transduction. For this reason, new investigations that allow the piezoelectric and electrokinetic behavior of the bone to be studied were undertaken [26].

A clinical study demonstrated that a local electromagnetic field accelerates the healing process after bone fracture [26]. Therefore, an article by Demiray and Dost [27] began new research concerning the effect of the electromagnetic field on interior injury to bone. In another article, Ramtani [26] presented a mathematical model relating to the benefit of the electric field in the reparair and maintenance of the solid matrix of bone. Furthermore, there are studies concerning the electrical behavior of bone tissue during the production of electric fields, and external electrical flow. Fukada and Yasuda [28] demonstrated that bone exhibits piezoelectric behavior, i.e. mechanical stress creates electric polarization (the indirect effect) and an external electric field causes strain (the converter effect). In addition, the properties of bones that produce piezoelectric potentials have been determined [29–33]. These data led to the development of mathematical models that include the effect of electromagnetic fields during the repair [34, 35] and remodeling [36] of bone. For example, Qu and Yu [34] developed a mathematical model (no spatial dimension) of the remodeling process and healing under the effect of mechanical loads and the use of electric charges. In this model, the higher the voltage applied to a bone after fracture, the lower the percentage of bone damage and micro damage in the few days after the stimulus. Similarly, during osteoporosis an electric field increases bone density over time. Huang et al. [37] established a hypothesis concerning the biological and biochemical pathways that activate cells, particularly osteocytes, during the imposition of an electric field. Furthermore, Qu and Yu [38] proposed a mathematical model that included mechanical loads and electromagnetic effects during the process of bone remodeling.

To date, there have been no phenomenological models concerning bone remodeling that have been tested and compared with purely mechanical models. Therefore, this article proposes a new electro-mechanical model relating to bone remodeling. To test its performance, various numerical experiments were carried out and compared with previous mechanical models. The electric model constants were obtained from relevant literature and numerical experimentation. From these assumptions it was concluded that the electric field can affect the distribution of mass, which originates from the remodeling process, under mechanical effects only. Using the remodeling model of Nackenhorst [39] as a starting point, it was demonstrated that the electric field can increase bone density and accelerate the process of apposition. Therefore, the model proposed herein can be used as a basis for further work concerning electrical effects in the maintenance of bone.

Methods

The electro-mechanical model

The electro-mechanical model of bone remodeling that involves mechanical and electrical stimuli can be written, hypothetically, as follows (1):

$$\frac{d\rho }{dt}=g_{\mathit{mech}}\left(\rho ,W_{\mathit{mec}}\left(\rho \right)\right)+g_{\mathit{elect}}\left(\rho ,W_{\mathit{elect}}\left(\epsilon (\rho ,f)\right)\right)$$

(1)

Where $g_{\mathit{mech}}\left(\rho ,W_{\mathit{mec}}\left(\rho \right)\right)$ is the well known mechanical stimulus described by Weinans [39], which depends on tissue density ( $\rho (x,y,z,t)$), and the work carried out by mechanical stress ( $W_{\mathit{mec}}\left(\rho \right)$) and $g_{\mathit{elect}}\left(\rho ,W_{\mathit{elect}}\left(\epsilon (\rho ,f)\right)\right)$ is the electrical stimulus that depends on density, frequency and the work carried out by the electric field ( $W_{\mathit{elect}}\left(\epsilon \left(\rho \right)\right)$). During this first approach, we consider that the two stimuli are added to determine the bone remodeling process. However, we will develop each of the terms that determine the electro-mechanical model throughout this article.

The mechanical model

Following the remodeling mechanical model described in [4, 39], the density variation over time depends on the mechanical stimulus that exists at every spatial point of the bone, which can be written as in [39] (2):

$$g_{\mathit{mech}}\left(\rho ,W_{\mathit{mec}}\left(\rho \right)\right)=k_1\left(\frac{W\left(\rho \right)}{W_{re{f}_m}}-1\right)$$

(2)

Where $\rho$is the bone density at each point in space ( $\rho (x,y,z,t)$), $W\left(\rho \right)$is the energy-strain per unit of volume due to mechanical stress, $k_1$ is a constant and $W_{\mathit{ref}}{}_m$ is the strain-energy (per unit of volume) of reference that sets the threshold for which they will perform deposition ( $W\left(\rho \right)/W_{re{f}_m}>1$) or absorption ( $W\left(\rho \right)/W_{re{f}_m}<1$) of the tissue [39], in the presence of mechanical stress. It should be noted that the energy-strain depends on the density and is given by (3):

$$W\left(\rho \right)=\frac{1}{2\rho }{\mathbf{\epsilon }}^{\mathbf{T}}\mathbf{C}\left(\rho \right)\mathbf{\epsilon }$$

(3)

Where $\mathbf{\epsilon }$ is the strain, in Voigt notation, of the strain tensor given in (4):

$$\mathbf{\epsilon }={\left[\begin{array}{ccc}\hfill {\epsilon }_{11}\hfill & \hfill {\epsilon }_{22}\hfill & \hfill {\gamma }_{12}\hfill \end{array}\right]}^T$$

(4)

That is a function of the displacements given by (5):

$$\mathbf{u}={\left[\begin{array}{cc}\hfill u_1\hfill & \hfill u_2\hfill \end{array}\right]}^T$$

(5)

$$\mathbf{\epsilon }={\mathbf{B}}_1\mathbf{u}=\left[\begin{array}{cc}\hfill \frac{\partial }{\partial x}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{\partial }{\partial y}\hfill \\ \hfill \frac{\partial }{\partial y}\hfill & \hfill \frac{\partial }{\partial x}\hfill \end{array}\right]\left[\begin{array}{c}\hfill u_1\hfill \\ \hfill u_2\hfill \end{array}\right]$$

(6)

In addition, $\mathbf{C}\left(\rho \right)$ is the matrix of linear elasticity. The matrix $\mathbf{C}\left(\rho \right)$ contains the Poisson module, which is usually considered constant, and Young's modulus that depends on the density by expression [4] (7):

$$E\left(\rho \right)=A{\rho }^n$$

(7)

where $A$ is a constant and $n$ establishes a relationship of power density that has been uncovered through experimentation [39].

By manipulating equation (7) we can obtain a dimensionless form of the density. This is easier to use with the aim of determining the modulus of elasticity. Multiplying the right side of (7) by $\left({\rho }_0^n/{\rho }_0^n\right)$ produces (8):

$$E\left(\rho \right)=A{\rho }_0^n{\left(\frac{\rho }{{\rho }_0}\right)}^n=E_0{\left(\frac{\rho }{{\rho }_0}\right)}^n=E_0{\lambda }^n$$

(8)

Where $E_0=A{\rho }_0^n$ and $\lambda =\rho /{\rho }_0$ are the elasticity modulus and the dimensionless density ratio, respectively. Therefore, the linear elasticity matrix can be expressed as (9):

$$\mathbf{C}\left(\rho \right)={\lambda }^n{\mathbf{C}}_0$$

(9)

Where ${\mathbf{C}}_0$ is the matrix of linear elasticity with constant coefficients, which depends on $E_0$ and $\nu$ only, and is given in the case of plane stress by:

$${\mathbf{C}}_0=\frac{E_0}{1-{\nu }^2}\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \nu \hfill & \hfill 0\hfill \\ \hfill \nu \hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{1}{2}\left(1-\nu \right)\hfill \end{array}\right]$$

(10)

Thus, the strain energy per unit of volume (3) can be expressed as (11):

$$W\left(\rho \right)=\frac{1}{2\rho }{\lambda }^n{\mathbf{\epsilon }}^{\mathbf{T}}{\mathbf{C}}_0\mathbf{\epsilon }=\frac{{\lambda }^n}{\rho }\left(\frac{{\mathbf{\epsilon }}^{\mathbf{T}}{\mathbf{C}}_0\mathbf{\epsilon }}{2}\right)=\frac{{\lambda }^n}{\rho }{\overline{U}}_{\mathit{mec}}$$

(11)

Where ${\overline{U}}_{\mathit{mec}}$ is the strain energy at each instant of time, which is calculated with the initial constant of the remodeling [39] problem only. Replacing these equations in (2), and with some algebraic manipulation, we produce (12):

$$g_{\mathit{mech}}\left(\rho ,W_{\mathit{mec}}\left(\rho \right)\right)=k_1\left(\frac{{\rho }_0}{{\rho }_0}\frac{{\lambda }^n{\overline{U}}_{\mathit{mec}}}{\rho W_{re{f}_m}}-1\right)=k_1\left(\left(\frac{{\rho }_0}{\rho }\right)\frac{{\lambda }^n{\overline{U}}_{\mathit{mec}}}{{\rho }_0{W}_{re{f}_m}}-1\right)=k_1\left({\lambda }^{n-1}\frac{{\overline{U}}_{\mathit{mec}}}{{\rho }_0{W}_{re{f}_m}}-1\right)$$

(12)

where we can define $U_{re{f}_m}={\rho }_0{W}_{re{f}_m}$. Therefore, we produce the following equation for the density ratio (13):

$$g_{\mathit{mech}}=k_1\left({\lambda }^{n-1}\frac{{\overline{U}}_{\mathit{mec}}}{U_{re{f}_m}}-1\right)$$

(13)

The momentum equation that establishes the internal stresses of a body is given by [40] (14):

$$\nabla ·\mathbf{\sigma }+\mathbf{b}=\mathbf{0}$$

(14)

Where the stress is given by (15):

$$\mathbf{\sigma }=\mathbf{C}\mathbf{\left(}\rho \mathbf{\right)}\mathbf{\epsilon }={\lambda }^n{\mathbf{C}}_0\mathbf{\epsilon }$$

(15)

Electrical model

This article proposes the inclusion of a hypothetical electrical term that can determine, in part, the process of bone remodeling. Therefore, the contribution of this stimulus can be written as (16):

$$g_{\mathit{elect}}\left(\rho ,W_{\mathit{elect}}\left(\epsilon (\rho ,f)\right)\right)=k_2\frac{W_{\mathit{elect}}\left(\epsilon \right(\rho ,f\left)\right)}{W_{re{f}_e}}$$

(16)

Where $\in (\rho ,f)$ is the electrical permittivity of bone tissue and depends on the density, ( $\rho$) where frequency ( $f$), $W_{\mathit{elect}}(\in (\rho ,f\left)\right)$ is the electrical energy per unit of volume, $k_2$ is a constant and $W_{\mathit{ref}}{}_e$ is the electric energy (per unit volume) of reference.

It should be noted that the term electric energy depends on permittivity, which in turn depends on the density and frequency. This energy term is given by (17):

$$W_{\mathit{elect}}(\in (\rho ,f\left)\right)=\frac{1}{2\rho }\in (\rho ,f){\mathbf{E}(x,y,z,t)}^2$$

(17)

Where $\mathbf{E}(x,y,z,t)$ is the electric intensity (electric field). However, $\epsilon (\rho ,f)$, electrical permittivity, depends on density and frequency, and is given by (18):

$$\in (\rho ,f)={\in }_0{\in }_r(\rho ,f)$$

(18)

Where ${\in }_0$ is the permittivity in the free space and ${\in }_r(\rho ,f)$ is the relative permittivity, which in turn is given by (19):

$${\in }_r(\rho ,f)=B{\rho }^m+\delta \left(f\right)$$

(19)

where $B$ is a constant and $m$ establishes a relationship of power density that can be proved through experimentation, and will be developed in the following sections. Furthermore, $\delta \left(f\right)$ is a function of frequency at which the electric field is applied.

As with the mechanical case, we can manipulate equation (19) so that it is expressed in terms of relative density, multiplying the first term by to produce $\left({\rho }_0^n/{\rho }_0^n\right)$ (20):

$${\in }_r(\rho ,f)=B{\rho }_0^m{\lambda }^m+\delta \left(f\right)={\in }_{\rho }{\lambda }^m+\delta \left(f\right)$$

(20)

where ${\in }_{\rho }$ is a constant for the power model of relative permittivity. Therefore, substituting (20) in (17), we obtain (21):

$$\begin{array}{l}{W}_{\mathit{elect}}(\in (\rho ,f\left)\right)=\frac{1}{2\rho }\frac{{\rho }_0}{{\rho }_0}{\in }_0({\in }_{\rho }{\lambda }^m+\delta (f\left)\right){\mathbf{E}(x,y,z,t)}^2\\ =\frac{1}{2{\rho }_0\lambda }{\in }_0({\in }_{\rho }{\lambda }^m+\delta (f\left)\right){\mathbf{E}(x,y,z,t)}^2=\frac{1}{{\rho }_0\lambda }({\in }_{\rho }{\lambda }^m+\delta (f\left)\right){\overline{U}}_{\mathit{elect}}\end{array}$$

(21)

where ${\overline{U}}_{\mathit{elect}}=\frac{1}{2}{\in }_0{\mathbf{E}(x,y,z,t)}^2$. Substituting equation (21) in (16) produces (22):

$$g_{\mathit{elect}}\left(\rho ,W_{\mathit{elect}}\left(\in \left(\rho \right)\right)\right)=k_2\frac{\frac{1}{{\rho }_0\lambda }({\in }_{\rho }{\lambda }^m+\delta (f\left)\right){\overline{U}}_{\mathit{elect}}}{W_{re{f}_e}}=k_2\frac{({\in }_{\rho }{\lambda }^m+\delta (f\left)\right)}{\lambda }\frac{{\overline{U}}_{\mathit{elect}}}{U_{re{f}_e}}$$

(22)

We have chosen $U_{re{f}_e}={\rho }_0{W}_{re{f}_e}$ as the reference value of electric energy.

However, Gauss's law for an electric field, with no internal loads, is given by (23):

$$\nabla ·\mathbf{D}={\rho }_e$$

(23)

where $\mathbf{D}=\in \mathbf{E}$ and ${\rho }_e$ is the density of the electric energy. In turn, the electric field can be expressed in terms of a quantity called electric potential or voltage, given by (24):

$$\mathbf{E}=-\nabla \varphi$$

(24)

where $\varphi$ is the electric potential.In summary, equation (1) can be written as (25):

$$\frac{d\rho }{dt}=k_1\left({\lambda }^{n-1}\frac{{\overline{U}}_{\mathit{mec}}}{U_{re{f}_m}}-1\right)+k_2\frac{({\in }_{\rho }{\lambda }^m+\delta (f\left)\right)}{\lambda }\frac{{\overline{U}}_{\mathit{elect}}}{U_{re{f}_e}}$$

(25)

Again, using the dimensionless form, we have (26):

$$\frac{d\lambda }{dt}=\frac{k_1}{{\rho }_0}\left({\lambda }^{n-1}\frac{{\overline{U}}_{\mathit{mec}}}{U_{re{f}_m}}-1\right)+\frac{k_2}{{\rho }_0}\frac{({\in }_{\rho }{\lambda }^m+\delta (f\left)\right)}{\lambda }\frac{{\overline{U}}_{\mathit{elect}}}{U_{re{f}_e}}$$

(26)

$$\frac{d\lambda }{dt}=k_{\mathit{mec}}\left({\lambda }^{n-1}\frac{{\overline{U}}_{\mathit{mec}}}{U_{re{f}_m}}-1\right)+k_{\mathit{elect}}\frac{({\in }_{\rho }{\lambda }^m+\delta (f\left)\right)}{\lambda }\frac{{\overline{U}}_{\mathit{elect}}}{U_{re{f}_e}}$$

(27)

where $k_{\mathit{mec}}=k_1/{\rho }_0$and $k_{\mathit{elect}}=k_2/{\rho }_0$ are mechanical and electrical constants that define the conversion rate of bone remodeling, dependent on the mechanical stress and the electrical potential, respectively.

Solution by the finite element method

To solve equation (14) we take its differential form to its weak [40] form, to obtain (27):

$$\underset{\Omega }{\int }\mathbf{\delta }\mathbf{\epsilon }{\lambda }^n{\mathbf{C}}_{\mathbf{0}}\mathbf{\epsilon }d\Omega -\underset{\Omega }{\int }\mathbf{\delta }\mathbf{u}·\mathbf{b}d\Omega -\underset{\Gamma }{\int }\mathbf{\delta }\mathbf{u}·\mathbf{t}d\Gamma =0$$

(28)

Where $\mathbf{\delta }\mathbf{\epsilon }$ and $\mathbf{\delta }\mathbf{u}$ are weighting functions, and $\mathbf{t}=\mathbf{\sigma }·\mathbf{n}$ is the traction on the boundary $\Gamma$ that serves as a frontier to the domain $\Omega$ (see Figure 1).

Figure 1

Domain and boundary conditions. The acronym on the index indicates whether it is a condition for the electric (elect) or mechanical (mec) case. Note that there are two types of boundaries (not necessarily equal) for each equation, the electric field and displacement (for the mechanical stress).

To calculate the approximate solution by finite element discretization, the displacement field approximates through [40] (28):

$$u(\mathbf{x},t)={\mathbf{N}}^T\mathbf{\left(}\mathbf{x}\mathbf{\right)}\mathbf{a}\left(t\right)$$

(29)

Where ${\mathbf{N}}^T\mathbf{\left(}\mathbf{x}\mathbf{\right)}$ is a row vector containing the shape functions used to approximate the displacement $\mathbf{a}\left(t\right)$ [40] in the nodes. Using the Galerkin method, we approximate the weighting functions in the same way as the displacement field, producing at the elementary level equation [40] (29):

$$\left(\underset{{\Omega }^e}{\int }\mathbf{B}{\lambda }^n{\mathbf{C}}_{\mathbf{0}}{\mathbf{B}}^{\mathbf{T}}d{\Omega }^e\right)\mathbf{a}+\underset{{\Omega }^e}{\int }\mathbf{N}\mathbf{b}d{\Omega }^e+B.C.=0$$

(30)

where $\mathbf{B}$ is the derivative operator (discrete version of operator ${\mathbf{B}}_1$ of equation (6)) that converts the displacements into strains (see [40]). This system of equations is completed by applying the Neumann and Dirichlet conditions suitable for solving the elastic problem (see Figure 1).

Similarly, for the electric case we can use the finite element method. For this purpose we use the weak expression of equation (23) in terms of electrical potential to produce (30):

$$\underset{\Omega }{\int }\left(\nabla w\right)\left({\in }_{\rho }{\lambda }^m+\delta \left(f\right)\right)\left(\nabla \varphi \right)d\Omega +\underset{\Omega }{\int }w{\rho }_{e}d\Omega -\underset{\Gamma }{\int }w\left(\left({\in }_{\rho }{\lambda }^m+\delta \left(f\right)\right)\left(\nabla \varphi \right)·\mathbf{n}\right)d\Gamma =0$$

(31)

where $w$ is the weighting function and $\mathbf{n}$ is the external normal to the domain of the problem. Again, for calculation of the approximate solution, the scalar field of the electric potential approaches through [40] (31):

$$\varphi (\mathbf{x},t)={\mathbf{N}}^T\mathbf{\left(}\mathbf{x}\mathbf{\right)}\mathbf{v}\left(t\right)$$

(32)

Where $\mathbf{v}$ is the potential at each node of the element, and $\mathbf{N}\mathbf{\left(}\mathbf{x}\mathbf{\right)}$ are the shape functions for scalar [40] problems. Weighting functions are chosen in the same manner as shape functions, which is through the method of (Bubnov-) Galerkin standard. Therefore, at the elementary level we have the following equation (32):

$$\left[\underset{\Omega }{\int }{\mathbf{B}}_{\mathbf{e}}\left({\in }_{\rho }{\lambda }^m+\delta \left(f\right)\right){\mathbf{B}}_e^{T}d{\Omega }^e\right]\mathbf{v}+\underset{\Omega }{\int }{\mathbf{N}}^{\mathbf{T}}{\rho }_{e}d\Omega +B.C.=0$$

(33)

where

$${\mathbf{B}}_{\mathbf{e}}=\left[\begin{array}{cc}\hfill \frac{\partial N_1}{\partial x}\hfill & \hfill \frac{\partial N_1}{\partial y}\hfill \\ \hfill \frac{\partial N_2}{\partial x}\hfill & \hfill \frac{\partial N_2}{\partial y}\hfill \\ \hfill \frac{\partial N_{\mathit{nod}}}{\partial x}\hfill & \hfill \frac{\partial N_{\mathit{nod}}}{\partial y}\hfill \end{array}\right]$$

(34)

where nod is the number of nodes of each element.

Solution to the equation of relative density

This article solves the equation of dimensionless density (or relative) by integrating equation (26) by the method of Euler. For this purpose we define (34):

$$\frac{d\lambda }{dt}=f(\mathbf{x},t_k,\lambda ,f)=k_{\mathit{mec}}\left({\lambda }^{n-1}\frac{{\overline{U}}_{\mathit{mec}}}{U_{re{f}_m}}-1\right)+k_{\mathit{elect}}\frac{({\in }_{\rho }{\lambda }^m+\delta (f\left)\right)}{\lambda }\frac{{\overline{U}}_{\mathit{elect}}}{U_{re{f}_e}}$$

(35)

Therefore, the forward Euler method is defined as [41] (35):

$${\lambda }_{k+1}={\lambda }_k+\Delta tf(\mathbf{x},t_k,{\lambda }_k)$$

(36)

where $\Delta t=t_{k+1}-t_k$ is the integration time interval and $k$ refers to the evaluation of the variable $\lambda$ at a specific time, i.e.: ${\lambda }_k=\lambda (\mathbf{x},t_k)$. This method has been used extensively in the prediction of bone density through remodeling [42]. The forward Euler method is of the first order and has the disadvantage of being unstable for large time intervals [41].

Euler's method was implemented in FORTRAN and was coupled with the elastic and electrical problems. For implementation we used an approach based on element (with an elemental average) [39, 43].

Computational model

The computational example to be solved is a 0.1x0.1 m square plate with non-uniform load at the top and restrictions of movement vertically at the bottom as presented in Figure 2. This example has been extensively studied in numerical tests concerning bone remodeling algorithms [42–44]. For these values we built a mesh of 80x80 elements (see Figure 3) in the horizontal and vertical directions, respectively, and integrated these with a time step of $\Delta t=0.100$. The total simulation time was 200. The initial condition was $\lambda =1.0$ and the limits used in the algorithm of bone remodeling were: $0.0125\le \lambda \le 2.175$.

Figure 2

Boundary conditions relating to bone remodeling.

Figure 3

Mesh used in the example. (a) 80x80.

The constants for the mechanical energy

For the mechanical case we used an initial Young’s modulus $E_0$, a Poisson's modulus $\nu$ and an initial dimensionless density ${\lambda }_0$ with values 64 MPa, 0.3 and 1.00, respectively [44]. The dimensionless parameters of the density equation were obtained from Nackenhorst (1997) [39, 44] and were $n=2$, $k_{\mathit{mec}}=0.3125day{s}^{-1}$ and $U_{re{f}_m}=800Pa$.

The mesh was produced using bilinear quadrilateral elements and four points of Gauss integration [40].

The constants for electric energy

In the electric case, there are numerous articles [45–48] that determine the major electrical properties of bone. In particular, graphs of the electric permittivity in function of frequency and density are provided in [46], as presented in Figures 4 and 5. Figure 4 presents data concerning bone density as a function of the relative permittivity on the frequency of 50 KHz. The article [46] used 40 bone tissue samples, and used the measurement of the density as a function of the permittivity. From the data, the mathematical model proposed for permittivity as a function of tissue density is given by a power equation, expressed as (36):

$${\in }_{\left(\rho \right)}(\rho ,50KHz)=B{\rho }^m$$

(37)

Figure 4

Bone density as a function of the relative permittivity at a frequency of 50 KHz. Taken from [46].

Figure 5

Graphic representation of relative permittivity (dimensionless) as a function of frequency (in Hz) for various densities of bone tissue. FC: Femoral head; FMC: femoral medial condyle, FLC: femoral lateral condyle FTM: femoral greater trochanter. Taken from [46].

Where $m=1.5486$ and $B=1050.0$; this is a function of the density of bone tissue. With the aim of using dimensionless density, we multiplied and divided by $\left({\rho }_0^n/{\rho }_0^n\right)$ to produce (37):

$${\in }_{\left(\rho \right)}\left(\rho \right)=1050.0{\rho }_0^{1.5486}{\lambda }^{1.5486}={\in }_{\rho }{\lambda }^{1.5486}$$

(38)

Where ${\in }_{\rho }=1050.0{\rho }_0^{1.5486}$ and depends on the initial density that considers the computational model.

Figure 5 presents the relative permittivity in terms of the function of the frequency for the density of the femoral head of $\rho =0.3737g.c{m}^{-3}$. With a greater frequency, permittivity decreases. To complement equation (37) (33) we must add a term that allows us to calculate the difference between the relative permittivity at a frequency of 50 kHz and any other frequencies used in the model; this is (38):

$$\begin{array}{l}{\delta }_{\left(f\right)}\left(f\right)={\in }_{r\left(f\right)}\left(0.3737g.c{m}^{-3},f\right)-{\in }_{r\left(50Khz\right)}\left(0.3737g.c{m}^{-3},f=50Khz\right)\\ {\delta }_{\left(f\right)}\left(f\right)={\in }_{r\left(f\right)}\left(0.3737g.c{m}^{-3},f\right)-228.681\end{array}$$

(39)

where ${\in }_{r\left(50Khz\right)}\left(0.3737g.c{m}^{-3},f=50Khz\right)=228.681$. For its part, the function ${\in }_{r\left(f\right)}\left(0.3737g.c{m}^{-3},f\right)$ can be obtained from Figure 5, so we get:

$${log}_{10}\left[{\in }_{r\left(f\right)}\left(0.3737g.c{m}^{-3},f\right)\right]=0.1407{\left[{log}_{10}\left(f\right)\right]}^2-1.8936{log}_{10}\left(f\right)+8.1505$$

(40)

In summary, the equation for the electric permittivity depending on the density and frequency is given by (see Figure 6):

$${\in }_r(\rho ,f)={\in }_{\rho }{\lambda }^{1.5486}+{\in }_{r\left(f\right)}-228.681$$

(41)

Figure 6

Graphic reconstruction of relative permittivity as a function of frequency and density.

For the remaining constant changes we made variations of $k_{\mathit{elec}}$ and used $U_{re{f}_e}=800Pa$, as in the mechanical case.

Results

Figure 7 presents various examples that were run for a frequency f = 50 KHz, and with variation in the constant K _elect . For these examples, a potential difference at various contour lines was imposed, according to the domain that was established in Figure 2. In the first example, a voltage of 100 was placed on the right side and no voltage at the bottom (Figure 7a). In the second example, a voltage of 100 V was placed on the right side and no voltage on the left (Figure 7b). In the third and fourth examples, the voltage was placed on the top and bottom, respectively (Figures 7c and 7d). In the other contours we placed null Neumann conditions.

Figure 7

Examples developed in the article; ϕ indicates voltage.

In the first example, the mechanical condition demonstrated in Figure 2 coupled with an electric field as observed in Figure 7a were used. The results are presented in Figure 8.

Figure 8

Results for the first example: the electric boundary condition in Figure 7a. The figure demonstrates the results for various values of k_elect.

Figure 8 presents the results for the relative density value of λ using an approach based on element, with an average at an elementary level of the mechanical stimulus, electrical and density values. In all figures (Figures 8, 9, 10, 11 and 12) black represents the maximum density value (2,175) and white represents the lowest value (0.0125). It should be noted that with Δt = 0.1 the pattern obtained is similar for all values of k _elect . It is observed that near the area of imposition of the mechanical stress there is a similar density distribution to a chess board, and far from the loading area there is the formation of three well-defined columns reminiscent of the cortical bone (Figure 8a and 8b). In cases where the constant increases (k _elect  = 1.4x10⁷ and above) there is the appearance of a high density area in the lower right region, from where a new column starts in the top domain (Figs  8c, d, e and f). In this area, near the lower right structure, there is the distribution and space formation of a chess board. In addition, in the upper part there are empty areas that generate ramifications in each of the columns supporting the load. Note that the imposition of the electric field defines a new topological structure in the domain, as presented on the right-hand side of the simulation results. This new structure is an additional column, which is generated by the potential and a support area of higher density in the lower right corner.

Figure 9

Results for the second example: varying the frequency in equations (35) and (36).

Figure 10

Results for the first example: the electric boundary condition in Figure 7b. The figure presents results for various values of k _elect .

Figure 11

Results for the first example: the electric boundary condition in Figure 7c. The figure presents results for various values of k _elect.

Figure 12

Results for the first example: the electric boundary condition in Figure 7 d. The figure demonstrates the results for various values of k _elect.

Figure 9 presents the results for various values of K _elect . In this case, a voltage of 100 V was imposed on the right side and a null voltage on the left side. Null Neumann conditions were imposed on remaining contours (see Figure 7b). It is noted that there were changes in topology; columns 2 and 3, from left to right in Figures 9c and 9d, are thicker and closer to one another. In addition, upper branches of greater density were created and an additional branch that begins from the last column to the right. In the case of Figure 9e, the formation of a non-defined region of the "chessboard" type on the right side of the domain is presented.

Figures 10 and 11 present results for various values of k _elect . In these cases a voltage of 100 V was imposed at the top and bottom, respectively. Null Neumann conditions were imposed on the remaining contours. Note that in these cases the mass distribution increases as k _elect increases. In the case of Figure 10e, the density increases in the upper part of the domain, so that the chessboard becomes more continuous in the central part of the domain. In Figures 10e and 11e the columns that are formed in each simulation have higher bandwidths than previous cases.

In the second example there is a variation of the frequency, while the value of k _elect remains constant at 7.0x10⁷. With the voltage configuration of Figure 7d, this is with a lower voltage of 100 V. Note in this case that at low frequencies the density and the amount of tissue deposited are greater than at high values. Figure 9a demonstrates that the frequency generates a total deposit of tissue. In Figure 12b the formation of a high density topology and continuing formation is apparent. Figures 12c, 12d and 12e present results comparable to those observed in previous cases. However, in Figure 12c, the columns are wider than in Figures 12d and 12e.

Discussion and conclusions

In this article several numerical examples were developed concerning bone remodeling of the plate during mechanical stress, assuming the imposition of an electric field in the domain. To calculate the mechanical and electrical stimulus of remodeling, and the evolution of density, the elemental approach was utilized. This pioneering article includes the electrical effect, previously designed by Weinans et al. [4], in a model of bone remodeling.

Comparable with previous articles, the results of the study presented herein demonstrate in the mechanical case formation of the columnar zone (of high density) in the area remote from the load and formation of the trabecular zone (of low density) in the area close to the load [17]. The results are similar to those obtained by Weinans et al. [4], Fernandez et al. [49] and Chen et al. [17]. When applying an electric field there is an increase in bone density and an alteration in the topology of the distribution of mass in the domain. In general, there is greater bone mass apposition in the domain. Therefore, the columns developed by the mechanical stress increase in size due to the electric field. In addition, a greater number of columns and localized compact zones may be observed. In the formulation the effect of electrical frequency has been included to allow increased apposition of mass at low frequency to be observed.