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.
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 contains the displacements and the electric potential . Using effective coefficients (elastic coefficients , piezoelectric coefficients , and dielectric coefficients ) and average state values (stress Cauchy tensor , strain tensor , electrical displacements , and electric fields ), the constitutive equations for homogeneous material can be expressed as [12–14]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) becomewhere , , and are defined as follows:The corresponding piezoelectric boundary integral equation for the boundary can be expressed as follows [16, 17]:where and are the load and field point, respectively, and is the free term coefficient at . and are the generalized traction and displacement vectors. The fundamental solution for elastic anisotropic piezoelectricity is given by [18, 19]where is a tensor function defined as [16, 20]and the fundamental solution for isThe 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 , 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 depends on the bone density and is constant [21, 22]. It is assumed thatwhere the constant β= 3 and M = 3790 x109 m8/kg2 s2.
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 areThe maximum values for the shear moduli are = 5.71 MPa, = 7.11 MPa, and = 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:with all other coefficients equal to zero. The variation of the density is calculated by adapting the equation (Eq. (3)) reported in [22]: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 is the strain energy density and is the local density. The corresponding equation iswhere B, D, and k are constants obtained from experiments and their values are 1.0 (gcm-3)2 (MPa per unit time), 60.0 (gcm-3)2 MPa (unit time), and 0.004 Jg-1, respectively, is the local density, and Ua represents the strain energy density (SED) modified herein to include the electric component:In order to obtain the real apparent density, it is assumed that the function in (14) is bounded by which is the maximum density of bone (cortical bone) and is the minimum density of bone (spongy bone). Resorption will take place when 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 .
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/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: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.
(a) Vertebra dimensions
(b) Vertebra discretization
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.
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.
(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., =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.
(a) Top view of the vertebra
(b) 3D view of the vertebra
(a) Upper view
(b) 3D view
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/cm3.
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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
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.