Skip to main content

Numerical test concerning bone mass apposition under electrical and mechanical stimulus


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.


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 [14]. 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 [68] 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 [1113] 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 [2933]. 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.


The electro-mechanical model

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

d ρ d t = g mech ρ , W mec ρ + g elect ρ , W elect ε ( ρ , f )

Where g mech ρ , W mec ρ is the well known mechanical stimulus described by Weinans [39], which depends on tissue density ( ρ ( x , y , z , t ) ), and the work carried out by mechanical stress ( W mec ( ρ ) ) and g elect ρ , W elect ε ( ρ , f ) is the electrical stimulus that depends on density, frequency and the work carried out by the electric field ( W elect ε ( ρ ) ). 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 mech ρ , W mec ρ = k 1 W ( ρ ) W r e f m 1

Where ρis the bone density at each point in space ( ρ ( x , y , z , t ) ), W ( ρ ) is the energy-strain per unit of volume due to mechanical stress, k 1 is a constant and W ref m is the strain-energy (per unit of volume) of reference that sets the threshold for which they will perform deposition ( W ( ρ ) / W r e f m > 1 ) or absorption ( W ( ρ ) / W r e 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 ( ρ ) = 1 2 ρ ε T C ( ρ ) ε

Where ε is the strain, in Voigt notation, of the strain tensor given in (4):

ε = ε 11 ε 22 γ 12 T

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

u = u 1 u 2 T
ε = B 1 u = x 0 0 y y x u 1 u 2

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

E ( ρ ) = A ρ n

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 ρ 0 n / ρ 0 n produces (8):

E ( ρ ) = A ρ 0 n ρ ρ 0 n = E 0 ρ ρ 0 n = E 0 λ n

Where E 0 = A ρ 0 n and λ = ρ / ρ 0 are the elasticity modulus and the dimensionless density ratio, respectively. Therefore, the linear elasticity matrix can be expressed as (9):

C ( ρ ) = λ n C 0

Where C 0 is the matrix of linear elasticity with constant coefficients, which depends on E 0 and ν only, and is given in the case of plane stress by:

C 0 = E 0 1 ν 2 1 ν 0 ν 1 0 0 0 1 2 1 ν

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

W ( ρ ) = 1 2 ρ λ n ε T C 0 ε = λ n ρ ε T C 0 ε 2 = λ n ρ U ¯ mec

Where U ¯ 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 mech ρ , W mec ρ = k 1 ρ 0 ρ 0 λ n U ¯ mec ρ W r e f m 1 = k 1 ρ 0 ρ λ n U ¯ mec ρ 0 W r e f m 1 = k 1 λ n 1 U ¯ mec ρ 0 W r e f m 1

where we can define U r e f m = ρ 0 W r e f m . Therefore, we produce the following equation for the density ratio (13):

g mech = k 1 λ n 1 U ¯ mec U r e f m 1

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

· σ + b = 0

Where the stress is given by (15):

σ = C ( ρ ) ε = λ n C 0 ε

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 elect ρ , W elect ε ( ρ , f ) = k 2 W elect ( ε ( ρ , f ) ) W r e f e

Where ( ρ , f ) is the electrical permittivity of bone tissue and depends on the density, ( ρ) where frequency ( f), W elect ( ( ρ , f ) ) is the electrical energy per unit of volume, k 2 is a constant and W 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 elect ( ( ρ , f ) ) = 1 2 ρ ( ρ , f ) E ( x , y , z , t ) 2

Where E ( x , y , z , t ) is the electric intensity (electric field). However, ε ( ρ , f ) , electrical permittivity, depends on density and frequency, and is given by (18):

( ρ , f ) = 0 r ( ρ , f )

Where 0 is the permittivity in the free space and r ( ρ , f ) is the relative permittivity, which in turn is given by (19):

r ( ρ , f ) = B ρ m + δ ( f )

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, δ ( f ) 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 ρ 0 n / ρ 0 n (20):

r ( ρ , f ) = B ρ 0 m λ m + δ ( f ) = ρ λ m + δ ( f )

where ρ is a constant for the power model of relative permittivity. Therefore, substituting (20) in (17), we obtain (21):

W elect ( ( ρ , f ) ) = 1 2 ρ ρ 0 ρ 0 0 ( ρ λ m + δ ( f ) ) E ( x , y , z , t ) 2 = 1 2 ρ 0 λ 0 ( ρ λ m + δ ( f ) ) E ( x , y , z , t ) 2 = 1 ρ 0 λ ( ρ λ m + δ ( f ) ) U ¯ elect

where U ¯ elect = 1 2 0 E ( x , y , z , t ) 2 . Substituting equation (21) in (16) produces (22):

g elect ρ , W elect ( ρ ) = k 2 1 ρ 0 λ ( ρ λ m + δ ( f ) ) U ¯ elect W r e f e = k 2 ( ρ λ m + δ ( f ) ) λ U ¯ elect U r e f e

We have chosen U r e f e = ρ 0 W r e 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):

· D = ρ e

where D = E and ρ 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):

E = ϕ

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

d ρ d t = k 1 λ n 1 U ¯ mec U r e f m 1 + k 2 ( ρ λ m + δ ( f ) ) λ U ¯ elect U r e f e

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

d λ d t = k 1 ρ 0 λ n 1 U ¯ mec U r e f m 1 + k 2 ρ 0 ( ρ λ m + δ ( f ) ) λ U ¯ elect U r e f e
d λ d t = k mec λ n 1 U ¯ mec U r e f m 1 + k elect ( ρ λ m + δ ( f ) ) λ U ¯ elect U r e f e

where k mec = k 1 / ρ 0 and k elect = k 2 / ρ 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):

Ω δ ε λ n C 0 ε d Ω Ω δ u · b d Ω Γ δ u · t d Γ = 0

Where δ ε and δ u are weighting functions, and t = σ · n is the traction on the boundary Γ that serves as a frontier to the domain Ω (see Figure 1).

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 ( x , t ) = N T ( x ) a ( t )

Where N T ( x ) is a row vector containing the shape functions used to approximate the displacement a ( t ) [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):

Ω e B λ n C 0 B T d Ω e a + Ω e N b d Ω e + B . C . = 0

where B is the derivative operator (discrete version of operator 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):

Ω w ρ λ m + δ ( f ) ϕ d Ω + Ω w ρ e d Ω Γ w ρ λ m + δ ( f ) ϕ · n d Γ = 0

where w is the weighting function and 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):

ϕ ( x , t ) = N T ( x ) v ( t )

Where v is the potential at each node of the element, and N ( x ) 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):

Ω B e ρ λ m + δ ( f ) B e T d Ω e v + Ω N T ρ e d Ω + B . C . = 0


B e = N 1 x N 1 y N 2 x N 2 y N nod x N nod y

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):

d λ d t = f ( x , t k , λ , f ) = k mec λ n 1 U ¯ mec U r e f m 1 + k elect ( ρ λ m + δ ( f ) ) λ U ¯ elect U r e f e

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

λ k + 1 = λ k + Δ t f ( x , t k , λ k )

where Δ t = t k + 1 t k is the integration time interval and k refers to the evaluation of the variable λ at a specific time, i.e.: λ k = λ ( 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 [4244]. 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 Δ t = 0.100 . The total simulation time was 200. The initial condition was λ = 1.0 and the limits used in the algorithm of bone remodeling were: 0.0125 λ 2.175 .

Figure 2
figure 2

Boundary conditions relating to bone remodeling.

Figure 3
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 ν and an initial dimensionless density λ 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 mec = 0.3125 d a y s 1 and U r e f m = 800 P a .

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 [4548] 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):

ρ ( ρ , 50 K H z ) = B ρ m
Figure 4
figure 4

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

Figure 5
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 ρ 0 n / ρ 0 n to produce (37):

ρ ( ρ ) = 1050.0 ρ 0 1.5486 λ 1.5486 = ρ λ 1.5486

Where ρ = 1050.0 ρ 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 ρ = 0.3737 g . 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):

δ f ( f ) = r ( f ) 0.3737 g . c m 3 , f r ( 50 K h z ) 0.3737 g . c m 3 , f = 50 K h z δ f ( f ) = r ( f ) 0.3737 g . c m 3 , f 228.681

where r ( 50 K h z ) 0.3737 g . c m 3 , f = 50 K h z = 228.681 . For its part, the function r ( f ) 0.3737 g . c m 3 , f can be obtained from Figure 5, so we get:

log 10 r ( f ) 0.3737 g . c m 3 , f = 0.1407 log 10 ( f ) 2 1.8936 log 10 ( f ) + 8.1505

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

r ( ρ , f ) = ρ λ 1.5486 + r ( f ) 228.681
Figure 6
figure 6

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

For the remaining constant changes we made variations of k elec and used U r e f e = 800 P a , as in the mechanical case.


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
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
figure 8

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

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.4x107 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
figure 9

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

Figure 10
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
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
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.0x107. 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.


  1. Ganong WF, William F: Fisiologia Médica. 2006, Editorial El Manual Moderno, Mexico DF, 373-386.

    Google Scholar 

  2. Cowin SC: Bone mechanics handbook. 2001, CRC Press, Washington

    Google Scholar 

  3. Jacobs CR, Simo JC, Beaupre GS, Carter DR: Adaptive bone remodeling incorporating simultaneous density and anisotropy considerations. J Biomech. 1997, 30: 603-613. 10.1016/S0021-9290(96)00189-3.

    Article  CAS  PubMed  Google Scholar 

  4. Weinans H, Huiskes R, Grootenboer H: The behavior of adaptive bone-remodeling simulation models. J Biomech. 1992, 25: 1425-1441. 10.1016/0021-9290(92)90056-7.

    Article  CAS  PubMed  Google Scholar 

  5. Wolff J: Das Gesetz der Transformation der Knochen. 1892, Hirschwald, Berlin, 1-139.

    Google Scholar 

  6. Frost HM: The laws of bone structure. 1964, CC Thomas, Springfield, IL

    Google Scholar 

  7. Frost HM: Mathematical elements of lamellar bone remodeling. 1964, Thomas, Springfield, IL

    Google Scholar 

  8. Frost HM: Vital biomechanics: proposed general concepts for skeletal adaptations to mechanical usage. Calcif Tissue Int. 1988, 42: 145-156. 10.1007/BF02556327.

    Article  CAS  PubMed  Google Scholar 

  9. Pauwels F: Gesammelte Abhandlungen Zur Funktionellen anatomie Des Bewegungsapparates. 1965, Springer-Verlag, Berlin, 543-549.

    Book  Google Scholar 

  10. Kummer B: Biomechanics of bone: mechanical properties, functional structure, functional adaptation. Biomechanics: its Foundations and Objectives. Edited by: Fung YC, Perrone N, Anliker M. 1972, Prentice-Hall, Englewood Cliffs, 237-272.

    Google Scholar 

  11. Cowin S: Wolffís law of trabecular architecture at remodeling equilibrium. J Biomech Eng. 1986, 108: 83-10.1115/1.3138584.

    Article  CAS  PubMed  Google Scholar 

  12. Cowin S, Hegedus D: Bone remodeling I: theory of adaptive elasticity. J Elast. 1976, 6: 313-326. 10.1007/BF00041724.

    Article  Google Scholar 

  13. Cowin S, Nachlinger RR: Bone remodeling III: uniqueness and stability in adaptive elasticity theory. J Elast. 1978, 8: 285-295. 10.1007/BF00130467.

    Article  Google Scholar 

  14. Hegedus D, Cowin S: Bone remodeling II: small strain adaptive elasticity. J Elast. 1976, 6: 337-352. 10.1007/BF00040896.

    Article  Google Scholar 

  15. Gupta S, New AMR, Taylor M: Bone remodelling inside a cemented resurfaced femoral head. Clin Biomech. 2006, 21: 594-602. 10.1016/j.clinbiomech.2006.01.010.

    Article  Google Scholar 

  16. Stülpner M, Reddy B, Starke G, Spirakis A: A three-dimensional finite analysis of adaptive remodelling in the proximal femur. J Biomech. 1997, 30: 1063-1066. 10.1016/S0021-9290(97)00074-2.

    Article  PubMed  Google Scholar 

  17. Jonkers I: Relation between subject-specific hip joint loading, stress distribution in the proximal femur and bone mineral density changes after total hip replacement. J Biomech. 2008, 41: 3405-3413. 10.1016/j.jbiomech.2008.09.011.

    Article  PubMed  Google Scholar 

  18. García J, Doblaré M, Cegoñino J: Bone remodelling simulation: a tool for implant design. Comput Mater Sci. 2002, 25: 100-114. 10.1016/S0927-0256(02)00254-9.

    Article  Google Scholar 

  19. Lian Z: Effect of bone to implant contact percentage on bone remodelling surrounding a dental implant. Int J Oral Maxillofacial Surgery. 2010, 39: 690-698. 10.1016/j.ijom.2010.03.020.

    Article  CAS  Google Scholar 

  20. Martin R: Targeted bone remodeling involves BMU steering as well as activation. Bone. 2007, 40: 1574-1580. 10.1016/j.bone.2007.02.023.

    Article  CAS  PubMed  Google Scholar 

  21. Hernandez C, Hazelwood S, Martin R: The relationship between basic multicellular unit activation and origination in cancellous bone. Bone. 1999, 25: 585-587. 10.1016/S8756-3282(99)00201-X.

    Article  CAS  PubMed  Google Scholar 

  22. Taylor D, Tilmans A: Stress intensity variations in bone microcracks during the repair process. J Theor Biol. 2004, 229: 169-177. 10.1016/j.jtbi.2004.03.014.

    Article  PubMed  Google Scholar 

  23. Hernandez C, Beaupre G, Carter D: A theoretical analysis of the changes in basic multicellular unit activity at menopause. Bone. 2003, 32: 357-363. 10.1016/S8756-3282(03)00037-1.

    Article  CAS  PubMed  Google Scholar 

  24. Peterson MC, Riggs MM: A physiologically based mathematical model of integrated calcium homeostasis and bone remodeling. Bone. 2010, 46: 49-63. 10.1016/j.bone.2009.08.053.

    Article  CAS  PubMed  Google Scholar 

  25. Buenzli PR, Pivonka P, Smith DW: Spatio-temporal structure of cell distribution in cortical Bone Multicellular Units: a mathematical model. Bone. 2010, 48 (4): 918-926.

    Article  PubMed  Google Scholar 

  26. Ramtani S: Electro-mechanics of bone remodelling. Int J Eng Sci. 2008, 46: 1173-1182. 10.1016/j.ijengsci.2008.06.001.

    Article  Google Scholar 

  27. Demiray H, Dost S: The effect of quadrupole on bone remodeling. Int J Eng Sci. 1996, 34: 257-268. 10.1016/0020-7225(95)00118-2.

    Article  Google Scholar 

  28. Fukada E, Yasuda I: On the piezoelectric effect of bone. J. Phys. Soc. Japan. 1957, 12: 1158-1162. 10.1143/JPSJ.12.1158.

    Article  Google Scholar 

  29. Aschero G, Gizdulich P, Mango F, Romano S: Converse piezoelectric effect detected in fresh cow femur bone* 1. J Biomech. 1996, 29: 1169-1174. 10.1016/0021-9290(96)00011-5.

    Article  CAS  PubMed  Google Scholar 

  30. Beck BR, Qin Y, Rubin C, McLeod K, Otter M: The Relationship of Streaming Potential Magnitude To Strain and Periosteal Modeling 567. Med & Sci Sports & Exercise. 1997, 29: 98-

    Article  Google Scholar 

  31. Gross D, Williams WS: Streaming potential and the electromechanical response of physiologically-moist bone. J Biomech. 1982, 15: 277-295. 10.1016/0021-9290(82)90174-9.

    Article  CAS  PubMed  Google Scholar 

  32. Hung C, Allen F, Pollack S, Brighton C: What is the role of the convective current density in the real-time calcium response of cultured bone cells to fluid flow?. J Biomech. 1996, 29: 1403-1409. 10.1016/0021-9290(96)84535-0.

    Article  CAS  PubMed  Google Scholar 

  33. Johnson MW, Chakkalakal DA, Harper RA, Katz JL: Comparison of the electromechanical effects in wet and dry bone. J Biomech. 1980, 13: 437-442. 10.1016/0021-9290(80)90037-8.

    Article  Google Scholar 

  34. Qu CY, Yu SW: The damage and healing of bone in the disuse state under mechanical and electro-magnetic loadings. Procedia Engineering. 2011, 10: 171-176.

    Article  Google Scholar 

  35. WANG E, ZHAO M: Regulation of tissue repair and regeneration by electric fields. Chinese Journal of Traumatology (English Edition). 2010, 13: 55-61.

    Google Scholar 

  36. Demiray H: Electro-mechanical remodelling of bones. Int J Eng Sci. 1983, 21: 1117-1126. 10.1016/0020-7225(83)90051-4.

    Article  Google Scholar 

  37. Huang CP, Chen XM, Chen ZQ: Osteocyte: the impresario in the electrical stimulation for bone fracture healing. Medical hypotheses. 2008, 70: 287-290. 10.1016/j.mehy.2007.05.044.

    Article  PubMed  Google Scholar 

  38. Qu C, Qin QH, Kang Y: A hypothetical mechanism of bone remodeling and modeling under electromagnetic loads. Biomaterials. 2006, 27: 4050-4057. 10.1016/j.biomaterials.2006.03.015.

    Article  CAS  PubMed  Google Scholar 

  39. Nackenhorst U: Numerical simulation of stress stimulated bone remodeling. Tech Mech. 1997, 17: 31-40.

    Google Scholar 

  40. Oñate E: Structural analysis with the finite element method. Linear statics. (Springer Verlag, 2009)

  41. Hoffman JD: Numerical methods for engineers and scientists. 1992, Mc Graw-Hill, New York

    Google Scholar 

  42. Beaupre G, Orr T, Carter D: An approach for time dependent bone modeling and remodelingótheoretical development. J Orthop Res. 1990, 8: 651-661. 10.1002/jor.1100080506.

    Article  CAS  PubMed  Google Scholar 

  43. Jacobs CR, Levenston ME, Beaupré GS, Simo JC, Carter DR: Numerical instabilities in bone remodeling simulations: the advantages of a node-based finite element approach. J Biomech. 1995, 28: 449-459. 10.1016/0021-9290(94)00087-K.

    Article  CAS  PubMed  Google Scholar 

  44. Chen G, Pettet G, Pearcy M, McElwain D: Comparison of two numerical approaches for bone remodelling. Medical engineering & physics. 2007, 29: 134-139. 10.1016/j.medengphy.2005.12.008.

    Article  CAS  Google Scholar 

  45. Sierpowska J: Prediction of mechanical properties of human trabecular bone by electrical measurements. Physiol Meas. 2005, 26: S119-10.1088/0967-3334/26/2/012.

    Article  CAS  PubMed  Google Scholar 

  46. Sierpowska J: Electrical and dielectric properties of bovine trabecular boneórelationships with mechanical properties and mineral density. Phys Med Biol. 2003, 48: 775-10.1088/0031-9155/48/6/306.

    Article  CAS  PubMed  Google Scholar 

  47. Gabriel C, Gabriel S, Corthout E: The dielectric properties of biological tissues: I. Literature survey. Phys Med Biol. 1996, 41: 2231-10.1088/0031-9155/41/11/001.

    Article  CAS  PubMed  Google Scholar 

  48. Sierpowska J: Interrelationships between electrical properties and microstructure of human trabecular bone. Phys Med Biol. 2006, 51: 5289-10.1088/0031-9155/51/20/014.

    Article  CAS  PubMed  Google Scholar 

  49. Fernández J, García-Aznar J, Martínez R, Viaño J: Numerical analysis of a strain-adaptive bone remodelling problem. Comput Methods Applied Mech Eng. 2010, 199: 1549-1557. 10.1016/j.cma.2010.01.005.

    Article  Google Scholar 

Download references


This work was financially supported by Division de Investigación de Bogotá, of Universidad Nacional de Colombia, under the modality: Fortalecimiento a grupos de investigación y creación artística con proyección nacional. ALIANZAS DE GRUPOS. The Project title is “Programa para la obtención y validación de parámetros numéricos utilizados en modelos matemáticos y simulación de sistemas y procesos biológicos mediante el diseño y caracterización de modelos celulares in vitro”, Hermes code 12873.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Diego A Garzón-Alvarado.

Additional information

Competing interests

The authors have no competing interests.

Authors&amp;#x2019; contributions

Work concerning production of the manuscript, modelling and numerical simulation was carried out by each of the authors. All authors read and approved the final manuscript.

Authors’ original submitted files for images

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Garzón-Alvarado, D.A., Ramírez-Martínez, A.M. & Cardozo de Martínez, C.A. Numerical test concerning bone mass apposition under electrical and mechanical stimulus. Theor Biol Med Model 9, 14 (2012).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: