A Simplified Scheme for Piezoelectric Anisotropic Analysis in Human Vertebrae Using Integral Methods

Abstract

This paper outlines a computational model for the analysis of the piezoelectric behaviour of the vertebral body remodelling process. Particular attention is paid to the algorithms for the simulation of the stress energy density for each point of the geometry and the distribution of the density in the bone. In addition, the model takes into account the piezoelectric effect and the anisotropy (transversal isotropy) of the bone. A model for internal anisotropic piezoelectric bone remodelling of a human vertebra is discussed in detail. The model consists of the implementation of an algorithm which includes the elastic and electric variables in a single equation using boundary element method. The presented results show a good agreement with biological data and the model does not include any electric additional charge.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Biomechanics is the application of mechanics to biology, trying to characterize the behaviour of living tissues and organs from a structural point of view. It also studies the changes due to different aspects and proposes methods for the modelling of the tensions in biological materials [1]. Living tissues are the most interesting known materials in terms of their structure and properties. The internal and external bone structure can be modified by inducing an electric potential as a result of loads acting on the bones. The relation between the formation and deformation of bones and the tension that generates electric potential makes these bones considered mechanically and electrically deformable, heterogeneous, and anisotropic [2, 3].

All properties as a whole are complex to study and model. The main purpose of the electric response of the bone during the application of loads or stresses is to show how the potential could influence the response to physiological functional forms [4]. Hence, even if the electromechanical interaction has been experimentally studied on different materials [5–7], the mathematical model, in biological terms, presented by many authors generally does not consider both the piezoelectric effect and transversal isotropic behaviour of the bone at the same time [8].

In this paper, it is assumed that bone electric properties influence the internal remodelling process of bones. How the density of bones grows and changes its form by the mechanical and electrical stimulus as an anisotropic material is explained. A three-dimensional (3D) mathematical model that introduces both mechanical and electric stimulus into the same equation is proposed. The formulation is based on the boundary element method (BEM) for homogeneous and anisotropic material. The piezoelectric 3D model includes the combination of the electric and mechanical properties and responses in an efficient computational algorithm. Its application to the internal bone remodelling process is discussed in detail and examples are presented.

2. Materials and Methods

2.1. Boundary Element Formulation

Piezoelectricity is commonly described as the constitutive coupling of electrical and mechanical fields (Coulomb’s law) [11]. The generalized displacement vector $U_k$ contains the displacements $u_k$ and the electric potential $\varphi$. Using effective coefficients (elastic coefficients $C_{ijkl}$, piezoelectric coefficients $e_{ikl}$, and dielectric coefficients ${\epsilon }_{ij}$) and average state values (stress Cauchy tensor ${\sigma }_{ij}$, strain tensor ${\gamma }_{kl}$, electrical displacements $D_i$, and electric fields $E_l$), the constitutive equations for homogeneous material can be expressed as [12–14]$$\begin{array}{}\text{(1)}& {\sigma }_{ij}=C_{ijkl}{\gamma }_{kl}-e_{lij}{E}_l,\text{(2)}& D_i=e_{ikl}{\gamma }_{kl}-{\epsilon }_{il}{E}_l\end{array}$$For piezoelectric materials, the problem can be formulated as in the elastic case using an extended displacement and traction vector with the electric potential and the electric charge [15]. The elastic and electric variables are combined in a single constitutive equation by introducing uppercase subscripts (J, K) which range from 1 to 4 and lowercase subscripts (i, l) which range from 1 to 3. Consequently, (1)-(2) become$$\begin{array}{}\text{(3)}& {\Sigma }_{iJ}=C_{iJKl}{Z}_{Kl}\end{array}$$where ${\Sigma }_{iJ}$, $Z_{Kl}$, and $C_{iJKl}$ are defined as follows:$$\begin{array}{c}\text{(4)}& {\Sigma }_{iJ}=\left\{\begin{array}{cc}{\sigma }_{ij},\hfill & j=J=\mathrm{1,2},\mathrm{3}\hfill \\ D_i,\hfill & J=\mathrm{4}\hfill \end{array}\right.Z_{Kl}=\left\{\begin{array}{cc}{\gamma }_{kl},\hfill & k=K=\mathrm{1,2},\mathrm{3}\hfill \\ -E_l,\hfill & K=\mathrm{4}\hfill \end{array}\right.\\ \text{(5)}& C_{iJKl}≔\left\{\begin{array}{cc}{C}_{ijkl},\hfill & j=J=\mathrm{1,2},3,\hspace{0.17em}\hspace{0.17em}K=\mathrm{1,2},3\hfill \\ e_{lij},\hfill & j=J=\mathrm{1,2},3,\hspace{0.17em}\hspace{0.17em}K=4\hfill \\ e_{ikl},\hfill & J=4,\hspace{0.17em}\hspace{0.17em}k=K=\mathrm{1,2},3\hfill \\ -{\epsilon }_{il},\hfill & J=4,\hspace{0.17em}\hspace{0.17em}K=4\hfill \end{array}\right.\end{array}$$The corresponding piezoelectric boundary integral equation for the boundary $\Gamma$ can be expressed as follows [16, 17]:$$\begin{array}{}\text{(6)}& C_{KJ}\left(\xi \right)U_J\left(\xi \right)={{\displaystyle \int }}_{\Gamma }^{}{U}_{KJ}^{\ast }\left(x,\xi \right)T_J\left(x\right)d\Gamma -{{\displaystyle \int }}_{\Gamma }^{}{T}_{KJ}^{\ast }\left(x,\xi \right)U_J\left(x\right)d\Gamma ,\end{array}$$where $\xi$ and $x$ are the load and field point, respectively, and $C_{KJ}$ is the free term coefficient at $\xi \in \Gamma$. $T_J$ and $U_J$ are the generalized traction and displacement vectors. The fundamental solution for elastic anisotropic piezoelectricity is given by [18, 19]$$\begin{array}{}\text{(7)}& U_{MK}^{\ast }=\frac{\mathrm{1}}{\mathrm{8}{\pi }^{\mathrm{2}}r}{{\displaystyle \oint }}_{\mathrm{0}}^{\mathrm{2}\pi }{\left(M_{MK}^{zz}\right)}^{-\mathrm{1}}\left(z\left(\varphi \right)\right)d\varphi \end{array}$$where $M_{MK}^{ab}$ is a tensor function defined as [16, 20]$$\begin{array}{}\text{(8)}& M_{MK}^{ab}=C_{iMKl}{a}_i{b}_l\end{array}$$and the fundamental solution for ${\mathrm{T}}_{\mathrm{M}\mathrm{J}}$ is$$\begin{array}{}\text{(9)}& T_{MJ}^{\ast }=C_{iJKl}{U}_{MK,l}^{\ast }{n}_i\end{array}$$The implementation of the anisotropic piezoelectric solution using the BEM requires the calculation of the fundamental solutions of the problem. The fundamental solutions for this specific problem cannot generally be solved analytically. Hence, the Radon Transform and the Dual Reciprocity Method (DRM) are combined in a framework which includes elastic and electric effects. This allows effectively solving the problem by applying the classical equations of elasticity extended to include the electric variables. It should be noted that the discretization and the equations of the BEM are restricted to the boundary of the problem. No volume discretization is required using this approach. For more details on the formulation and its implementation the reader is referred to [17].

2.2. Piezoelectric Anisotropic 3D Model for Bone Remodelling

The model for piezoelectric bone remodelling is based on the following physical observations of the process:

(i)

The bone is a piezoelectric material that is associated with the presence of oriented fibrous proteins such as collagen. These show hexagonal symmetry behaviour.

(ii)

The piezoelectric material depends on the anisotropy of the material, which means that the elastic matrix should be considered in the model.

(iii)

The internal and external bone structure can be modified by an induced electric potential which results from loads acting on the bones. A corresponding relationship exists between the formation and deformation of bone and the stress generated potential (piezoelectric effect).

(iv)

Young’s modulus $E$, as well as piezoelectric and dielectric matrices, depends on the bone density, while Poisson’s ratio ν is constant.

From this initial assumption, the adaptability of bones can be modelled by the electric and mechanical reactions which change the apparent density. This is known as internal bone remodelling process. The constitutive tensor matrix is composed of an elastic matrix considering a transversal isotropic material, a piezoelectric matrix, and a dielectric matrix. In the proposed piezoelectric remodelling model, spongy and compact bone is considered. The spongy bone is defined as a homogeneous and isotropic material where $E$ depends on the bone density and $\nu$ is constant [21, 22]. It is assumed that$$\begin{array}{c}\text{(10)}& E=M{\rho }^{\beta }\nu =\mathrm{0.3}\end{array}$$where the constant β= 3 and M = 3790 x10⁹ m⁸/kg² s².

In the following, compact bone is modelled as an anisotropic material where the constitutive tensor matrix is aligned with the distribution of the osteons in the compact bone [23]. The material parameters are$$\begin{array}{c}\text{(11)}& E_{\mathrm{1}}=E_{\mathrm{2}}=\mathrm{2314}{\rho }^{\mathrm{1.57}}{E}_{\mathrm{3}}=\mathrm{2065}{\rho }^{\mathrm{3.09}}\\ G_{\mathrm{12}}=\frac{G_{\mathrm{12}max}{\rho }^{\mathrm{2}}}{{\rho }_{max}^{\mathrm{2}}}\\ G_{\mathrm{23}}=\frac{G_{\mathrm{23}max}{\rho }^{\mathrm{2}}}{{\rho }_{max}^{\mathrm{2}}}\\ G_{\mathrm{31}}=\frac{G_{\mathrm{31}max}{\rho }^{\mathrm{2}}}{{\rho }_{max}^{\mathrm{2}}}\\ {\nu }_{\mathrm{12}}=\mathrm{0.4}\\ {\nu }_{\mathrm{23}}=\mathrm{0.25}\\ {\nu }_{\mathrm{31}}=\mathrm{0.25}\end{array}$$The maximum values for the shear moduli $G$ are $G_{\mathit{\text{12}}max}$ = 5.71 MPa, $G_{\mathit{\text{23}}max}$ = 7.11 MPa, and $G_{\mathit{\text{31}}max}$ = 6.58 MPa.

Transversal isotropic behaviour is assumed, i.e., E₁ = E₂ = E_p, E₃ = E_t, ν₁₂ = ν₂₁ = ν_p, ν₃₁ = ν₃₂ = ν_tp, and ν₁₃ = ν₂₃ = ν_pt. Hence, the coefficients of the constitutive matrix can be simplified into the following form:$$\begin{array}{c}\text{(12)}& C_{\mathrm{11}}={\mathrm{C}}_{\mathrm{22}}={\mathrm{E}}_p\left(\mathrm{1}-{\nu }_{pt}{\nu }_{tp}\right)\zeta {\mathrm{C}}_{\mathrm{33}}={\mathrm{E}}_t\left(\mathrm{1}-{\nu }_p^{\mathrm{2}}\right)\zeta \\ {\mathrm{C}}_{\mathrm{12}}={\mathrm{E}}_p\left({\nu }_p+{\nu }_{pt}{\nu }_{tp}\right)\zeta \\ {\mathrm{C}}_{\mathrm{13}}={\mathrm{E}}_p\left({\nu }_{tp}+{\nu }_p{\nu }_{tp}\right)\zeta ={\mathrm{E}}_t\left({\nu }_{pt}-{\nu }_p{\nu }_{pt}\right)\zeta \\ C_{\mathrm{44}}=\mathrm{2}{G}_{\mathrm{12}}\\ C_{\mathrm{55}}=\mathrm{2}{G}_{\mathrm{31}}\\ C_{\mathrm{66}}=\mathrm{2}{G}_{\mathrm{23}}\\ \zeta =\frac{\mathrm{1}}{\mathrm{1}-{\nu }_p^{\mathrm{2}}-{\mathrm{2}\nu }_{pt}{\nu }_{tp}-\mathrm{2}{\nu }_p{\nu }_{pt}{\nu }_{tp}}\\ \mathrm{C}\mathrm{21}=\mathrm{C}\mathrm{12}\\ \mathrm{C}\mathrm{31}=\mathrm{C}\mathrm{23}=\mathrm{C}\mathrm{32}=\mathrm{C}\mathrm{13}\end{array}$$with all other coefficients equal to zero. The variation of the density is calculated by adapting the equation (Eq. (3)) reported in [22]:$$\begin{array}{}\text{(13)}& {\Sigma }_{iJ}={\left(\frac{\rho }{{\rho }^{\ast }}\right)}^{\beta }{C}_{iJKl}{Z}_{Kl}.\end{array}$$where ${\rho }^{\ast }$ is equal to 1 and the variation of $\rho$ is a first-order equation [24]. In the remodelling theory developed in [25], the change in the bone density ($\rho$) is expressed as a function of the mechanical stimulus[21], where $U(\sigma ,ϵ\left(u\right))$ is the strain energy density and $\rho$ is the local density. The corresponding equation is$$\begin{array}{}\text{(14)}& \frac{d\rho }{dt}=B\left(\frac{U_a}{\rho }-k\right)-D{\left(\frac{U_a}{\rho }-k\right)}^{\mathrm{2}}\hspace{1em}{\rho }_{tb}<\rho \le {\rho }_{cb}\end{array}$$where B, D, and k are constants obtained from experiments and their values are 1.0 (gcm^-3)² (MPa per unit time), 60.0 (gcm^-3)² MPa ${}^{\text{-2}}$ (unit time)${}^{\text{-1}}$, and 0.004 Jg^-1, respectively, $\rho$ is the local density, and U_a represents the strain energy density (SED) modified herein to include the electric component:$$\begin{array}{}\text{(15)}& U_a\left(\sigma ,ϵ\left(u\right)\right)=\frac{\mathrm{1}}{\mathrm{2}}{Z}_{iJ}{\Sigma }_{iJ}\end{array}$$In order to obtain the real apparent density, it is assumed that the function in (14) is bounded by ${\rho }_{cb}$ which is the maximum density of bone (cortical bone) and ${\rho }_{tb}$ is the minimum density of bone (spongy bone). Resorption will take place when $d\rho /dt$ is negative; otherwise bone deposition occurs.

2.3. Algorithm for Piezoelectric Bone Remodelling

The formulation outlined in the previous sections was implemented into a BEM framework. Available BEM codes [20, 26] for static elasticity problems were adapted by the authors to include piezoelectricity and dynamics. The DRM is used to provide accurate results for the simulation of the dynamic piezoelectric bone remodelling process. The density is calculated at the nodes and projected on all nodes of the model in every time step. It was observed that the zones with more density are associated with bone deposition, while the bone loss is associated with reabsorption.

Figure 1 shows a flowchart of the proposed computational algorithm. The algorithm begins with the definition of the 3D geometry, i.e., the nodes and boundary elements, and the boundary conditions for load, electric charge, and density. The density is used to calculate the most important qualitative properties of anisotropic piezoelectricity and elastic, piezoelectric, and dielectric matrix. These properties are changed according to the density for each element in the geometry. The fundamental solutions are calculated numerically using the Radon Transformation [14, 17, 20]. Finally, the density can be computed from the stress energy density at every point in the geometry and all values in the system can be updated, until the end of the simulation time ${\mathrm{T}}_{\mathrm{s}\mathrm{i}\mathrm{m}}$.

[figure omitted; refer to PDF]

2.4. Application

In the following, the application of the presented BEM approach to piezoelectric bone remodelling is shown. The model corresponds to a 3D simulation of the vertebra (lumbar region) and the visualization of the numerical results was carried out using the GiD software[27]. The biomechanical model is presented in [28]. The authors analyzed a two-dimensional model with loads from the daily activity and represented the initial model conditions of the bone remodelling postsurgery after an implant. The example presented in the following considers a vertebra without apophysis. The bone is simulated with an elastic isotropic matrix when the density is between 0 and 1.08 g/cm³ (spongy bone (10)) and with a hexagonal symmetry of a crystal (transversal isotropic) for densities >1.08 g/cm³ (compact bone (11)). The values of the piezoelectric and dielectric constants are given below:$$\begin{array}{c}\text{(16)}& e=\left[\begin{array}{cccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{17.88215}& \mathrm{3.57643}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{3.57643}& -\mathrm{17.88215}& \mathrm{0}\\ \mathrm{1.050765}& \mathrm{1.050765}& \mathrm{1.87209}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\times {\mathrm{10}}^{\mathrm{9}}C/{mm}^{\mathrm{2}}\epsilon =\left[\begin{array}{ccc}88.54& 0& 0\\ 0& 88.54& 0\\ 0& 0& 106.248\end{array}\right]\times {10}^{-12}{C}^2/{Nmm}^2\end{array}$$The values of the piezoelectric and dielectric matrices change according to the density in each element on the mesh and following the form presented in [22, 24]. The current approach also incorporates the anisotropic effect on the remodelling equation (see (13)). The discretization of the boundary of the vertebra with quadrilateral boundary elements is shown in Figure 2. A total of 494 quadrilateral boundary elements are used to represent the boundary of the vertebra.

[figures omitted; refer to PDF]

Flexion and extension are common in lumbar spine. Lateral flexion is free at the atlantooccipital joint. Rotation is smaller at the lumbar region. Hence the boundary conditions shown in Figure 3 are assumed, where u, v, and w correspond to the directions of the x-, y-, and z-axis.

[figure omitted; refer to PDF]

The Neumann conditions of the model were taken from [28, 29] and shown in Figure 4. In the young state, the shape of the load distribution is a symmetrically concave parabola with ps_max=4.8 N/mm2 (two sides) and ps_min=1.6 N/mm2 (center) on the top and pi_max =4.525 N/mm2 (two sides) and pi_min =1.325 N/mm2 (center) at the bottom.

[figures omitted; refer to PDF]

3. Results and Discussion

The results obtained with the current approach are compared to results presented in the literature where the bone density was determined from a tomographic image.

The simulation ended after 85 days of simulation time in one CPU i5 8GB RAM (i.e., ${\mathrm{T}}_{\mathrm{s}\mathrm{i}\mathrm{m}}$=85 days) and the initial density of the bone was 0.8 g/cm³. Figure 5(a) shows a cortical bone density in the range of 1.41 to 1.7379 g/cm³ on the left and right walls which represent cortical bone. On the walls to the central upper and lower part of the vertebra the range is between 0.27 and 1.11g/cm³, corresponding to spongy or trabecular bone of normal vertebral bodies. From Figure 6 it can also be noticed that the density is higher at the vertebra boundary and that minor ossification occurs in the middle of the vertebral body. The simulation results in Figure 5 show a very good agreement between the observations of the biological images and the simulations. It is interesting to note that the simulated bone remodelling gives values of cortical and trabecular bone which are close to the one that can be observed at the dry bone data. However, it is difficult to make an objective comparison.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

Figure 7 shows the predicted density distribution along the vertical axis (sagittal plane). The cortical bone can clearly be seen in the periphery whereas the trabecular bone can be seen in the middle of the vertebral body. These observations are in good agreement with those made from the medical tomographic data shown in Figure 8. However, Figure 5(b) does not show a strong apophysis. In this area, the bone is trabecular with a density of 0.27 to 0.92g/cm³.

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

Although the numerical simulation does not include an additional electric charge surface, it has been shown that the algorithm incorporating the electric variables yields the expected results which are in good agreement with the typical medical data. Therefore, the algorithm could be used to simulate the addition of electrical charges and develop a pacing protocol that can reduce the medical problems associated with bone loss.

4. Conclusions

A BEM piezoelectric model based on a number of assumptions and a density of strain energy equation with the addition of the electric conditions has been developed.

The internal remodelling algorithm described in this paper details the equation and the dynamical processes which are used in order to simulate the incremental changes in density. The methodology is valid on a vertebral body, with similar characteristics based on the physical assumptions of the daily activity. The piezoelectric bone remodelling model of the vertebra accurately predicts the anisotropy in the cortical bone. The mechanical and electric stimuli agree with how the bone changes, even though an electric surface of the model was not used.

The numerical model presented could be used to study the influence of electric stimulation on osteoporosis. Future studies should analyze the total anisotropic remodelling and the addition of the surface charge on the vertebral body. These can bring new findings that contribute to the development of new therapies or implants which can help in keeping the bone’s integral shape and decrease bone loss and fracture risks.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The first author would like to acknowledge the financial support of the Venezuelan Academy of Natural Sciences and ALFA II ELBENET (Europe Latin America Boundary Element Network) for supporting her doctoral studies and CIMNE Classrooms for their cooperation and support for GiD. The authors extend their gratitude to Professor Gernot Beer for giving permission to use the BEM code and for his invaluable advices during the research.

References

[1] D. Schneck, J. Bronzino, Biomechanics Principles and Applications,DOI: 10.1201/9781420040029, 2003.

[2] A. Gjelsvik, "Bone remodeling and piezoelectricity - I," Journal of Biomechanics, vol. 6 no. 1, pp. 69-77, DOI: 10.1016/0021-9290(73)90039-0, 1973.

[3] A. Marino, R. Becker, "Origen of the piezoelectric effect in bone," Calcified Tissue Research, pp. 177-180, 1971.

[4] L. A. Norton, "Effects of variable frequency on in vitro bone growth in an electric field," Journal of Dental Research, vol. 53 no. 3, pp. 771-771, DOI: 10.1177/00220345740530035101, 1974.

[5] R. O. Becker, "The bioelectric factors in amphibian-limb regeneration," The Journal of Bone & Joint Surgery—American Volume, vol. 43, pp. 643-656, DOI: 10.2106/00004623-196143050-00002, 1961.

[6] M. H. Shamos, "The origin of bioelectric effects in mineralized tissues," Journal of Dental Research, vol. 44 no. 6, pp. 1114-1115, DOI: 10.1177/00220345650440060901, 1965.

[7] E. Fukada, I. Yasuda, "On the piezoelectric effect of bone," Journal of the Physical Society of Japan, vol. 12 no. 10, pp. 1158-1162, DOI: 10.1143/jpsj.12.1158, 1957.

[8] M. Gupta, D. K. Sinha, "On a mathematical model in Biopiezoelectricity Indian," J. pure appl. Math, pp. 11-522, 1980.

[9] E. Perilli, I. Parkinson, K. Reynolds, "Micro-CT examination of human bone: from biopsies towards the entire organ," Ann. Ist. Super. Sanità, 2012.

[10] H. Visser, N. Motsitsi, "Percutaneous Trans Ilio-Sacral Biopsy of the Spine in sacral tuberculosis," Internet Journal of Orthopedic Surgery, vol. 8, 2007.

[11] M. Kögl, M. L. Bucalem, "A family of piezoelectric MITC plate elements," Computers & Structures, vol. 83 no. 15-16, pp. 1277-1297, DOI: 10.1016/j.compstruc.2004.04.025, 2005.

[12] M. Denda, C.-Y. Wang, "3D {BEM} for the general piezoelectric solids," Computer Methods Applied Mechanics and Engineering, vol. 198 no. 37-40, pp. 2950-2963, DOI: 10.1016/j.cma.2009.04.014, 2009.

[13] Y. González, V. Duarte, F. Nieto, M. Cerrolaza, "The boundary element method in computational biomechanics: some applications," Issues on Mechanical and Civil Engineering, 2012.

[14] V. Duarte, Y. González, M. Cerrolaza, "Boundary element simulation of bone tissue," International Journal of Biomedical Engineering and Technology, vol. 5 no. 2-3, pp. 211-228, DOI: 10.1504/IJBET.2011.039199, 2011.

[15] F. García-Sánchez, A. Sáez, J. Domínguez, "Anisotropic and piezoelectric materials fracture analysis by BEM," Computers & Structures, vol. 83 no. 10-11, pp. 804-820, DOI: 10.1016/j.compstruc.2004.09.010, 2005.

[16] L. Gaul, M. Kögl, M. Wagner, Boundary Element Methods for Engineers and Scientists,DOI: 10.1007/978-3-662-05136-8, 2003.

[17] M. Cerrolaza, V. Duarte, D. Garzón-Alvarado, "Analysis of Bone Remodeling Under Piezoelectricity Effects Using Boundary Elements," Journal of Bionic Engineering, vol. 14 no. 4, pp. 659-671, DOI: 10.1016/S1672-6529(16)60432-8, 2017.

[18] N. A. Schclar, Topics in Engineering. Anisotropic Analysis Using Boundary Elements, vol. 20, 1994.

[19] J. A. Sanz, M. P. Ariza, J. Dominguez, "Three-dimensional BEM for piezoelectric fracture analysis," Engineering Analysis with Boundary Elements, vol. 29 no. 6, pp. 586-596, DOI: 10.1016/j.enganabound.2004.12.014, 2005.

[20] K. Thoeni, Efiziente Berechnung anisotroper Fundamentallösungen für die Methode der Randelemente[Master's, thesis], 2005.

[21] H. Weinans, R. Huiskes, H. J. Grootenboer, "The behavior of adaptive bone-remodeling simulation models," Journal of Biomechanics, vol. 25 no. 12, pp. 1425-1441, DOI: 10.1016/0021-9290(92)90056-7, 1992.

[22] J. R. Fernández, J. M. García-Aznar, R. Martínez, "Piezoelectricity could predict sites of formation/resorption in bone remodelling and modelling," Journal of Theoretical Biology, vol. 292, pp. 65-81, DOI: 10.1016/j.jtbi.2011.09.032, 2012.

[23] V. Báča, D. Kachlík, Z. Horák, J. Stingl, "The course of osteons in the compact bone of the human proximal femur with clinical and biomechanical significance," Surgical and Radiologic Anatomy, vol. 29 no. 3, pp. 201-207, DOI: 10.1007/s00276-007-0192-6, 2007.

[24] R. Martínez, Numerical Analysis and Simulations in Bone Remodeling Models [Ph.D. thesis], 2010.

[25] J. Li, H. Li, L. Shi, A. S. L. Fok, C. Ucer, H. Devlin, K. Horner, N. Silikas, "A mathematical model for simulating the bone remodeling process under mechanical stimulus," Dental Materials, vol. 23 no. 9, pp. 1073-1078, DOI: 10.1016/j.dental.2006.10.004, 2007.

[26] G. Beer, I. Smith, C. Duenser, The Boundary Element Method with Programming, 2008.

[27] GiD, CIMNE, "GiD the personal pre and post processor," .

[28] J. D. Jovanovic, M. L. Jovanovic, Biomechanical model of vertebra based on bone remodeling, Facta Universitatis: Medicine and Biology, 2004.

[29] Z. Xinghua, G. He, Z. Dong, G. Bingzhao, "A study of the effect of non-linearities in the equation of bone remodeling," Journal of Biomechanics, vol. 35 no. 7, pp. 951-960, DOI: 10.1016/S0021-9290(02)00028-3, 2002.