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

Authors: Vannessa de J. Duarte; Klaus Thoeni; Diego Garzón-Alvarado; Miguel Cerrolaza
Journal: Mathematical Problems in Engineering
 (2018)

---

## 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.


---

## Body

## 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 Uk contains the displacements uk and the electric potential ϕ. Using effective coefficients (elastic coefficients Cijkl, piezoelectric coefficients eikl, and dielectric coefficients εij) and average state values (stress Cauchy tensor σij, strain tensor γkl, electrical displacements Di, and electric fields El), the constitutive equations for homogeneous material can be expressed as [12–14](1)σij=Cijklγkl-elijEl,(2)Di=eiklγkl-εilElFor 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(3)ΣiJ=CiJKlZKlwhere ΣiJ, ZKl, and CiJKl are defined as follows:(4)ΣiJ=σij,j=J=1,2,3Di,J=4ZKl=γkl,k=K=1,2,3-El,K=4(5)CiJKl≔Cijkl,j=J=1,2,3,K=1,2,3elij,j=J=1,2,3,K=4eikl,J=4,k=K=1,2,3-εil,J=4,K=4The corresponding piezoelectric boundary integral equation for the boundary Γ can be expressed as follows [16, 17]:(6)CKJξUJξ=∫ΓUKJ∗x,ξTJxdΓ-∫ΓTKJ∗x,ξUJxdΓ,where ξ and x are the load and field point, respectively, and CKJ is the free term coefficient at ξ∈Γ. TJ and UJ are the generalized traction and displacement vectors. The fundamental solution for elastic anisotropic piezoelectricity is given by [18, 19](7)UMK∗=18π2r∮02πMMKzz-1zϕdϕwhere MMKab is a tensor function defined as [16, 20](8)MMKab=CiMKlaibland the fundamental solution for TMJ is(9)TMJ∗=CiJKlUMK,l∗niThe 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 modulusE, 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 whereE depends on the bone density and ν is constant [21, 22]. It is assumed that(10)E=Mρβν=0.3where the constantβ= 3 andM = 3790 x109m8/kg2s2.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(11)E1=E2=2314ρ1.57E3=2065ρ3.09G12=G12maxρ2ρmax2G23=G23maxρ2ρmax2G31=G31maxρ2ρmax2ν12=0.4ν23=0.25ν31=0.25The maximum values for the shear moduli G are G12max= 5.71 MPa, G23max= 7.11 MPa, and G31max= 6.58 MPa.Transversal isotropic behaviour is assumed, i.e.,E1= E2= Ep, E3= Et,ν12=ν21=νp,ν31=ν32=νtp,andν13=ν23=νpt. Hence, the coefficients of the constitutive matrix can be simplified into the following form:(12)C11=C22=Ep1-νptνtpζC33=Et1-νp2ζC12=Epνp+νptνtpζC13=Epνtp+νpνtpζ=Etνpt-νpνptζC44=2G12C55=2G31C66=2G23ζ=11-νp2-2νptνtp-2νpνptνtpC21=C12C31=C23=C32=C13with all other coefficients equal to zero. The variation of the density is calculated by adapting the equation (Eq. (3)) reported in [22]:(13)ΣiJ=ρρ∗βCiJKlZKl.where ρ∗ is equal to 1 and the variation of ρ is a first-order equation [24]. In the remodelling theory developed in [25], the change in the bone density (ρ) is expressed as a function of the mechanical stimulus[21], where U(σ,ϵu) is the strain energy density and ρ is the local density. The corresponding equation is(14)dρdt=BUaρ-k-DUaρ-k2ρtb<ρ≤ρcbwhereB, D, andk are constants obtained from experiments and their values are1.0 (gcm-3)2 (MPa per unit time),60.0 (gcm-3)2MPa-2 (unit time)-1, and0.004 Jg-1, respectively, ρ is the local density, andUa represents the strain energy density (SED) modified herein to include the electric component:(15)Uaσ,ϵu=12ZiJΣiJIn order to obtain the real apparent density, it is assumed that the function in (14) is bounded by ρcb which is the maximum density of bone (cortical bone) and ρtb is the minimum density of bone (spongy bone). Resorption will take place when dρ/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.Figure1 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 Tsim.Figure 1
Flowchart of the proposed computational algorithm for BEM piezoelectric bone remodelling.


### 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 between0 and1.08 g/cm3 (spongy bone (10)) and with a hexagonal symmetry of a crystal (transversal isotropic) for densities >1.08 g/cm3 (compact bone (11)). The values of the piezoelectric and dielectric constants are given below:(16)e=00017.882153.5764300003.57643-17.8821501.0507651.0507651.87209000×109C/mm2ε=88.5400088.54000106.248×10-12C2/Nmm2The 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.Figure 2
Vertebra discretization (dimension in mm), mesh of 1160 nodes, and 494 quadrilaterals boundary elements.
(a)
Vertebra dimensions (b)
Vertebra discretizationFlexion 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 Figure3 are assumed, where u, v, and w correspond to the directions of thex-,y-, andz-axis.Figure 3
Boundary conditions of the vertebra.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 withps_max=4.8 N/mm2 (two sides) andps_min=1.6 N/mm2 (center) on the top andpi_max =4.525 N/mm2 (two sides) andpi_min =1.325 N/mm2 (center) at the bottom.Figure 4
Biological load on the vertebra.
(a)
Tridimensional view (b)
Bidimensional view


## 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 Uk contains the displacements uk and the electric potential ϕ. Using effective coefficients (elastic coefficients Cijkl, piezoelectric coefficients eikl, and dielectric coefficients εij) and average state values (stress Cauchy tensor σij, strain tensor γkl, electrical displacements Di, and electric fields El), the constitutive equations for homogeneous material can be expressed as [12–14](1)σij=Cijklγkl-elijEl,(2)Di=eiklγkl-εilElFor 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(3)ΣiJ=CiJKlZKlwhere ΣiJ, ZKl, and CiJKl are defined as follows:(4)ΣiJ=σij,j=J=1,2,3Di,J=4ZKl=γkl,k=K=1,2,3-El,K=4(5)CiJKl≔Cijkl,j=J=1,2,3,K=1,2,3elij,j=J=1,2,3,K=4eikl,J=4,k=K=1,2,3-εil,J=4,K=4The corresponding piezoelectric boundary integral equation for the boundary Γ can be expressed as follows [16, 17]:(6)CKJξUJξ=∫ΓUKJ∗x,ξTJxdΓ-∫ΓTKJ∗x,ξUJxdΓ,where ξ and x are the load and field point, respectively, and CKJ is the free term coefficient at ξ∈Γ. TJ and UJ are the generalized traction and displacement vectors. The fundamental solution for elastic anisotropic piezoelectricity is given by [18, 19](7)UMK∗=18π2r∮02πMMKzz-1zϕdϕwhere MMKab is a tensor function defined as [16, 20](8)MMKab=CiMKlaibland the fundamental solution for TMJ is(9)TMJ∗=CiJKlUMK,l∗niThe 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 modulusE, 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 whereE depends on the bone density and ν is constant [21, 22]. It is assumed that(10)E=Mρβν=0.3where the constantβ= 3 andM = 3790 x109m8/kg2s2.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(11)E1=E2=2314ρ1.57E3=2065ρ3.09G12=G12maxρ2ρmax2G23=G23maxρ2ρmax2G31=G31maxρ2ρmax2ν12=0.4ν23=0.25ν31=0.25The maximum values for the shear moduli G are G12max= 5.71 MPa, G23max= 7.11 MPa, and G31max= 6.58 MPa.Transversal isotropic behaviour is assumed, i.e.,E1= E2= Ep, E3= Et,ν12=ν21=νp,ν31=ν32=νtp,andν13=ν23=νpt. Hence, the coefficients of the constitutive matrix can be simplified into the following form:(12)C11=C22=Ep1-νptνtpζC33=Et1-νp2ζC12=Epνp+νptνtpζC13=Epνtp+νpνtpζ=Etνpt-νpνptζC44=2G12C55=2G31C66=2G23ζ=11-νp2-2νptνtp-2νpνptνtpC21=C12C31=C23=C32=C13with all other coefficients equal to zero. The variation of the density is calculated by adapting the equation (Eq. (3)) reported in [22]:(13)ΣiJ=ρρ∗βCiJKlZKl.where ρ∗ is equal to 1 and the variation of ρ is a first-order equation [24]. In the remodelling theory developed in [25], the change in the bone density (ρ) is expressed as a function of the mechanical stimulus[21], where U(σ,ϵu) is the strain energy density and ρ is the local density. The corresponding equation is(14)dρdt=BUaρ-k-DUaρ-k2ρtb<ρ≤ρcbwhereB, D, andk are constants obtained from experiments and their values are1.0 (gcm-3)2 (MPa per unit time),60.0 (gcm-3)2MPa-2 (unit time)-1, and0.004 Jg-1, respectively, ρ is the local density, andUa represents the strain energy density (SED) modified herein to include the electric component:(15)Uaσ,ϵu=12ZiJΣiJIn order to obtain the real apparent density, it is assumed that the function in (14) is bounded by ρcb which is the maximum density of bone (cortical bone) and ρtb is the minimum density of bone (spongy bone). Resorption will take place when dρ/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.Figure1 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 Tsim.Figure 1
Flowchart of the proposed computational algorithm for BEM piezoelectric bone remodelling.


## 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 between0 and1.08 g/cm3 (spongy bone (10)) and with a hexagonal symmetry of a crystal (transversal isotropic) for densities >1.08 g/cm3 (compact bone (11)). The values of the piezoelectric and dielectric constants are given below:(16)e=00017.882153.5764300003.57643-17.8821501.0507651.0507651.87209000×109C/mm2ε=88.5400088.54000106.248×10-12C2/Nmm2The 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.Figure 2
Vertebra discretization (dimension in mm), mesh of 1160 nodes, and 494 quadrilaterals boundary elements.
(a)
Vertebra dimensions (b)
Vertebra discretizationFlexion 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 Figure3 are assumed, where u, v, and w correspond to the directions of thex-,y-, andz-axis.Figure 3
Boundary conditions of the vertebra.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 withps_max=4.8 N/mm2 (two sides) andps_min=1.6 N/mm2 (center) on the top andpi_max =4.525 N/mm2 (two sides) andpi_min =1.325 N/mm2 (center) at the bottom.Figure 4
Biological load on the vertebra.
(a)
Tridimensional view (b)
Bidimensional view


## 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.,Tsim=85 days) and the initial density of the bone was 0.8 g/cm3. Figure 5(a) shows a cortical bone density in the range of 1.41 to 1.7379 g/cm3 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/cm3, 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.Figure 5
Top view and 3D view of the density on the vertebral body (g/cm3).
(a)
Top view of the vertebra (b)
3D view of the vertebraFigure 6
Tomographic image of the vertebral body after [9].
(a)
Upper view (b)
3D viewFigure7 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/cm3.Figure 7
Distribution of the density along the vertical Z axis on the middle vertebral body (g/cm3).Figure 8
Saggital plane view of the vertebral body from [10].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.



---

*Source: 2783867-2018-06-19.xml*