Simulation of Subject Specific Bone Remodeling

1. Introduction

Vertebral compression fractures are very common low-bone-density related conditions among the elderly and have a significant adverse health effect [1]. According to a survey by the National Osteoporosis Foundation, osteoporosis related fractures cost an estimated $19 billion in 2005 and the cost is predicted to go up to $25.3 billion by the year of 2025 [2]. The fractures affects up to 25% of postmenopausal women [3]. Men are equally affected by this, but because of greater peak bone mass, it occurs about 10 years later than women [4]. The thoraco-lumbar region in the vertebra is the most common region where fractures occur, with T8, T12, L1 and L4 being the most commonly fractured segments [5]. The vertebral compression fracture can be avoided by preventative medication and physical therapy if the trend of the bone structure change is known. However, directly monitoring the density and structural change of patient’s bones is practically impossible based on current technology. Examination of the patients’ bones at regular intervals is costly, and not very helpful: the patients have to wait for months or years to know the effect of their lifestyles on their vertebra, while at this point any bone changes might already have occurred and damages made. Instead, simulation, as an cost-effective and non-intrusive method, is a better way to “know” the bone change in advance based on realistic input of loads and patients conditions, and thus prevent bone fractures occurring.

Changes in the cancellous bone structure have been studied since JuliusWolff first discovered the adaptive nature of the bone [6]. Different simulation techniques and models have been proposed [7, 8]. However, in the previous studies standard load which is within the nominal range experienced by the human vertebra were used. In reality this load varies for different human subjects for various activities performed. This was a limitation in terms of the ability to measure accurate loads. With technologies like the state-of-the-art motion capture systems combined with musculoskeletal modeling and simulation softwares, accurate loads on different segments of the body can be calculated. This helps in more accurate simulation results.

This study aims at developing a simulation tool to predict the long term effects on the micro-architecture of the vetebral cancellous bone, given the person’s present daily lifestyle, within a reasonable time frame (hours or days). Modules are created based on different purposes of simulation, and thesemodules are integrated into one efficient simulationmethod (Figure 1) that covers from the acquisition of individual motion data to the prediction of the changes under various physical loads. The process starts with scanning of the vertebral bodies to obtain μCT images. The μCT images aremicro-computed tomography (μCT) images generated using x-rays. The result is a 2D cross section of a 3D object which shows the internal structure of the 3D object. In this particular case, pre-scanned μCT images were obtained. These images are processed to separate the bone and the marrow. A 3D model of the bone is constructed from the processed images. Quantifying parameters are calculated on the 3D model of the bone. After calculating the loads on the vertebrae caused by different activities, finite element methods are used to calculate the stresses within the structure caused by those loads. This triggers the remodeling process which continues until a steady state is reached and there are no more significant changes in the structure. If the method can differentiate and quantify the vertebra bone structure changes under different types of physical loads, then the method can be used to produce desirable physical therapy against the bone fractures and density loss. The details of the method are discussed in the following sections.

Figure 1.

Schematic of proposed remodeling cycle

2. Material and methods

Bone by its nature is dynamic. It undergoes constant structural changes as a result of adaptation and remodeling. The purpose of remodeling is to prevent the accumulation of damage, adapt the internal architecture to external loads [9]. To effectively study the remodeling phenomenon, it is important to know the components of the bone and the remodeling process [10, 11].

2.1. Image processing

Tomography is widely used in medical radiology to produce noninvasive, diagnostic, cross-sectional images of bone and tissue structure in human patients. As the input to our modeling, we use a series of 2D μCT images, i.e., tomographic images whose resolution is in the micrometer range [12]. Figure 2(a) shows a sample image which is 16 bit gray-scale image, equivalent to histological sections used for microscopic slides. In order to effectively reconstruct a 3D geometry of the bone, it is important to identify the various regions in the gray scale image. The image shows bone, marrow and any artifact that might have been introduced due to the instrument and the imaging technique or post-processing technique used by the image reconstruction software. The goal of this section is to segment the bone region in the μCT images and generate a 3D geometry from these images.

Figure 2.

a) OriginalμCT scan (b) Binary image after thresholding the ROI

An important step in the analysis of μCT images is segmentation. This process separates the bone from the surrounding marrow substance. Because of the difference in the density between the bone and the marrow, the region occupied by the calcified bone is lighter compared to the region occupied by the marrow. An automated global thresholding method proposed by Otsu [13], which chooses the threshold to minimize the interclass variance of the black and white pixels in a grayscale image, was used to separate the bone from marrow. Since the focus of this project is on the remodeling and adaptation of the cancellous bone, a region of interest was chosen from the center of the vertebra. This removes the cortical shell and any artifacts present at the edge during scanning. The diameter of the ROI was arbitrarily chosen to be 480 pixels. To maintain the aspect ratio (diameter/height) of that of a single vertebral body (1.33), a total of 360 images (480/360=1.33) were chosen from the entire set. These cored images were then thresholded to separate the bone matrix form the surrounding marrow (Figure 2(b)).

Once all the images were thresholded and segmented, the next step is to construct the 3D geometry from them. 3D geometry can be represented in many ways. One way is to store a list of vertices on the surface and a list of indices that form connected triangular facets from those vertices. This is an efficient and the most commonly used file format to store 3D geometries. An isosurface algorithm is used to extract connected triangular facet information. The result of this is a 3D surface geometry representing the cored section of the cancellous bone (Figure 3).

Figure 3.

Surface model of the segmented core

Another useful way to represent the 3D geometry is using voxels. A voxel is a 3D pixel element, i.e., a cube of user specified dimensions. The resolution of the voxel depends on the resolution and the slice thickness of the μCT image. If the μCT images are stacked on top of each other in a sequence, the slice number in the stack defines the Zcoordinate. A right hand coordinate systemis used throughout the document unless otherwisementioned. The Xand Yvalue of each pixel representing the bone in an image slice along with the slice number of that particular image in the stack, defines the X, Yand Zcoordinate of each voxel. An example of the voxel model is illustrated in Figure 4. Although the surface model is a concise and efficient way of representing the 3D geometries, the voxel model is useful in calculating the various parameters that quantify the structure as explained in the following section.

2.2. 3D thinning algorithm

Thinning can be considered as a type of data compression where the important geometric features of a structure are retained without having to store a large amount of data to represent the actual geometry. Because of this, thinning is used in 3D object and feature recognition methods. Recognizing the features of a 3D geometry helps in accurately measuring the change in those features. Thus thinning serves as an intermediate step in calculating the direct 3D quantifying parameters of the bone. Thinning generates a mid-axis for a rod-like structure and a mid-plane for a plate-like structure. The process involves systematically shaving off the surface voxels without changing the topology of the structure until the skeleton of the structure is reached, at which point no more surface voxel can be removed. The skeletons computed by various thinning methods produce spurious curves and surface branches [15, 16]. These short branches do not add any useful information. Pruning the unwanted short branches from both the curves and the surfaces with user controllable parameters is an useful feature, especially for complex structures like cancellous bone. Ju et al. [17] proposed an efficient pruning algorithm, which handles both the curve and the surface branches so that the final skeleton is both topology preserving and shape depicting. A slight variation of the method proposed by Ju et al. has been successfully implemented in MATLAB.

The overview of the process in given below,

1. Find the bounday voxel

2. Find the critical points in the set of boundary voxels

3. A point is a critical point if its removal changes the topology of the structure or if it is a curve or a surface end point

4. Remove the non-critical points from the set of bondary points

5. Repeate steps 1-3 untill there are no more boundary points to remove

Figure 4.

Voxel model of the segmented core

The next step in thinning is pruning. Noisy and redundant surface and curve edges are removed in this step. This works on the same principal as morphological 2D image operators, erosion and dilation. These two operations are performed on the set of object points or voxels. Similar to themorphological operation on 2D images, opening has the effect of removing small irregular protrusions and produces a smooth edge. In the complete algorithm to calculate the skeleton of an 3D object, the thinning and pruning processes are run alternatively on the object points [11, 17].

2.2.1. Validation of thinning algorithm

The performance of thinning algorithm was tested on shapes with varying structural complexity in structure. It was first tested on a simple CAD geometry of regular grid of rods of known dimension and thickness. As expected the algorithm generated a clean one voxel thick mid-axis. This is shown in Figure 5, where the gray region represents the CAD structure and the blue line is the calculated mid-axis. This structure was rotated by 5 and 30 around the X-axis, and the skeleton was calculated on both the geometries. The result on these structures is discussed in [11].

Figure 5.

CAD model and its mid-axis overlapped on it.

Next, the algorithm was tested on a CAD model with plates in it and μCT scans of Duocel open cell aluminum foam with 20 pores per inch. In terms of the structure metal foam falls between the regular CAD file and the cancellous bone. The results and explanation in given in [11]. Finally a cancellous bone sample was used to calculate mid-axes; the results of which are shown in Figure 6. This information will be used in the next section (2.3) to calculate parameters quantifying the structure.

Figure 6.

a) Bone sample and (b) its skeleton overlapped on it.

2.2.2. Parallelization of thinning algorithm

The thinning algorithms are iterative in nature, i.e., voxels are traversed and removed either layer by layer or sequentially. As the size or the resolution of the geometry increases, the time taken to calculate the skeleton increases. For a complex structure like the cancellous bone within a vertebra, it becomes practically impossible to calculate the skeleton on a single machine within a reasonable time frame (minutes or even hours). To get an approximation of the computational time required to do such a calculation, Figure 7 shows a plot of time required to do such a calculation for different size geometries. If the calculations were being done on an entire vertebra, the time taken would be too long to be of any practical use.

Figure 7.

Time taken to calculate the skeleton.

One solution to this problem is to share the data across multiple processors and do the calculations in a distributed way. The concept of parallel computation is not new [18, 19]. Various techniques exists to perform parallel computations on a given data set. However, these algorithms are specific to an application and the literature on parallel computation of medial axes for the purpose of quantification is scarce. A simple and efficient method was developed and implemented using MATLAB’s Parallel Computing Toolbox (PCT) and Distributed Computing Server toolboxes [14]. The concept is to split the data into smaller fragments which is customizable in terms of size and dimensions, send it to different processor cores, perform 3D thinning on the individual segments with the knowledge of how the segments are shared, collect the resulting data and combine the fragments to generate the final skeleton of the original structure. The details of the implementation are discussed in [11].

2.3. Quantification

Cancellous bone undergoes constant change in its micro-architecture. This change is more prominent in bone disorders like osteoporosis. As a result, the bone mass changes. Bone mass alone is not sufficient to quantify these changes. Clinical studies have shown that the cancellous structure experiences significant changes which affect the mechanical properties of the bone changes [20, 21]. This change in the micro-architecture of cancellous bone leads to decrease in bone strength and eventually lead to vertebral compression fractures [22]. A prominent change of plate-like structures to rod-like structures was shown by Ding et al. [21] after the age of 60 years. Loss of entire struts or conversion of plates to struts due to perforation was shown by Parfitt et al. [23]. These previous data show that apart from the overall bone mass, cancellous bone goes through important structural changes.

Apart from the visual assessment of the structural changes, different structural parameters are used to characterize these changes [24, 25]. In the past, these parameters were studied by the examination of the 2D cross-sections of cancellous bone biopsies. The 3D morphometric parameters are then derived from 2D images using stereological methods assuming a fixed structual model [24]. Recent advances in μCT have made it possible to acquire these parameters directly from 3D μCT image [26, 27] without any underlying model assumptions. The calculated midaxes can also be used to decompose the bone into rods and plates [17]. Further, the orientation of the rods can also be calculated. By decomposing the cancellous bone into rods and plates, the parameters can be calculated for the rods and plates seperately. A comparison between the indirect calculated values assuming fixed model structure and the model independent direct calculation is made in [27]. These comparisons show that making assumptions on the type of cancellous bone structure leads to errors in the calculated values. Hence, a model-independent calculation of the parameters was used.

2.3.1. Parameter calculations

Based on the literature search the following set of parameters were chosen to characterize the bone sample [24, 27]. They are: bone volume fraction, which is a ratio of bone volume to the total volume (BV/TV); surface density, which is the ratio of bone surface area to the bone volume (BS/BV); trabecular thickness (Tb.Th); trabecular separation (Tb.Sp); and the trabecular number (Tb.N). The cancellous thickness, separation and number are measured directly from the 3D geometry with the method proposed by Hildebrand et al. [26] as described later in the section.

The bone volume (BV) is the volume enclosed by the bone surface. This does not include the marrow space. Total volume (TV) is the volume of the bone including the marrow space. Since the bone geometry is irregular, it is difficult to calculate the volume. To overcome this problem, the bone geometry is split into smaller known geometries (a tetrahedron in this case) and add the volume of each element [11]. Tetrahedral mesh generation algorithm [28] was used to divide the bone in to smaller elements. The bone surface area is calculated by adding the area of all the triangular facets generated by the isosurface function [10, 11].

Before calculating the other three parameters, the skeleton of the structure is calculated as described in section 2.2. The thickness calculation uses a volume-based approach rather than the conventional surface based approach [26]. According to the volume-based approach, the local thickness (τ) is defined as below. If Vis the set of object points, then the local thickness measured at point p∈ Vis the diameter of the largest sphere centered at point qcontaining the point pand completely inside the structure at the same time (Figure 8). The mean thickness of the structure is defined as the arithmetic mean of all the local thicknesses. To implement this on a discrete 3D volume, the Euclidean distance from a point (q) to the nearest background point is assigned to that point. This is calculated using a fast distancemap algorithm proposed by Maurer et al. [29]. This distance is equivalent to the radius rof the sphere centered at that point q. The local thickness measured at point qis given by Equation 1 where, sph(q, r) is the set of points within the sphere centered at qwith a radius r.

$$\tau (q)=2\ast \max (r|sph(q,r)\subseteq V)$$E1

Figure 8 illustrates this concept. The blue line is the object boundary, red line is the mid-axis of the object. If pis a point within the structure and p∈ sph(q, r) then the local thickness of point pis given by 2 ∗ r, the diameter of the largest sphere which encloses point p. Since the points on the mid-axes are the center of largest diameter sphere that can be fit within the structure, the local thickness of all the points on the mix-axes are calculated. The mean thickness of the structure is the average of all the local thicknesses calculated for that structure.

Figure 8.

Thickness calculation of a structure.

Calculating the trabecular separation (Tb.Sp) is similar except the maximum diameter sphere is fit within the marrow space. For trabecular number (Tb.N), the sphere fit within the marrow space is extended all the way to the mid-axes of the adjacent trabeculae.

2.3.2. Algorithm verification

The quantifying algorithm explained and implemented in the previous sections was verified on three different types of structures, a 3D CAD model of regular grid and known dimension with only rods, a similar CAD model with both rods and plates, and finally on μCT images of open cell metal foam structure with 20 pores per inch. The 3D CAD models have an aspect ratio of 1.3 which is the same as vertebral bone [30]. Since the parameters are measured on digitized images, the resolution of the μCT images have a significant effect on the calculated values. The 3D CAD models were imported into MATLAB with varying resolutions, and the set of parameters were calculated for each resolution. The error vs. resolution graph is plotted for each of the parameters for the 3D CAD model (Figure 9). As can be seen from these graphs, as the resolution increases (in other words, the size of voxels decreases) the error becomes small. The other cause of error due to digitization is how well the voxels are aligned with the surface boundary of the 3D geometry.

Figure 9.

Parameters of CAD model in its original orientation (a) Trabecular thickness (b) Trabecular separation (c) Trabecular number (d) Bone volume fraction (e) Bone surface density.

The quantifying parameters were measured on this CAD after rotating it around the x-axis by 5°and then by 30°. The same procedure was repeated on a CAD model with both rods and plates. The results of this is discussed in [11].

The third structure used for verification was a metal foam sample scanned at a resolution of 28.7μm/pixel. The average calculated values of the parameters for this structure are shown in Table 1.

Parameters	Calculated values	Measured values	% Error
rTb.Th (mm)	0.3276	0.3198	-2.4405%
pTb.Th (mm)	0.4618	0.4558	-1.3163%
Tb.Sp (mm)	0.7118	0.7312	2.6532%
Tb.N (mm−1 )	1.0999	0.9991	-10.08%
BV/TV (%)	0.1301	0.1319	1.3646%
BS/BV (mm−1 )	8.3297	8.4213	1.0877%

Table 1.

Calculated parameters of the metal foam sample.

To find the accurateness of the results, struts were randomly picked from the metal foam sample. The thickness of these struts were measured by hand by fitting spheres within the struts. To measure the angle, cylinders were fit to the strut and the angle of the axis of the cylinder was taken as the true angle of the strut. These values were then compared to the calculated values from the quantifying algorithm. The metal foam sample has 300 struts in it. To assess the results with a 95% confidence level and within a 5% confidence interval, the sample size was determined to be 169 [42]. Lillefors test for normality was performed on the data sets using the MATLAB’s Statistical toolbox which proved that the both of the data sets come from data which is normally distributed. Thus the verification process shows that the calculated values are within an acceptable range of ±5% of error. It was also concluded that the resolution and the orientation of the scans of the object have significant effect on the calculated values. As the resolution increases the error in the calculated values are generally reduced.

2.4. Lumbar load calculation

The skeleton is a load bearing structure. Skeletal segments are connected to each other by joints of various degrees of freedom. Each joint in the body experiences different amounts of load depending on the type of activity done by the human subject. This subject specific loading data can be used to predict the state of the bone structure with techniques that have been studied, tested and verified to comply to the rules of bone remodeling (2.5). The realistic loading data can be collected by surgically placing the load sensors into the subject’s body at a desired location [31, 32]. But this is not a practical solution to be used as a clinical tool. An alternative solution is to use motion data of the subject to drive a human body model in a musculoskeletal modeling and simulation software like lifemod [43].

2.4.1. Motion capture (MOCAP)

The first step is to capture the motion of a subject. The Vicon motion capture system was used for this purpose. The setup consists of 10 ViconMX cameras. A full body Plug-in-Gait marker set with 39 markers were used. 9mm retro-reflective markers were fixed on to the black jumpsuit, worn by the user, with velcro. The details of the marker set are given in [11]. For the preliminary study, to show the proof of concept, two different types of motion were chosen: walking, and rope jumping exercise. Basic subject specific data was recorded for each subject. Two subjects were used in this test. Once the subject were markered up, s/he was asked to perform the two tasks and the data was collected using Vicon Nexus software. Figure 10 shows a screen shot of the Vicon Nexus software illustrating the camera setup and the subject within the capture volume. A custom MATLAB script was written to convert the motion capture data from Vicon system into LifeMOD compatible format.

Figure 10.

Camera setup for motion capture.

2.4.2. Detailed lumbar model in LifeMOD

LifeMOD is a virtual human modeling and simulation software. It is built on MD ADAMS software which is a standard for mechanical system simulation. Given the basic information such as age, height, weight and gender, LifeMOD can generate a full body model with the basic bones, muscles and joints based on its anthropometric database. This basic model can be refined according to the users needs by adding more muscles and/or bone segments. It gives the user the ability to develop a full-body or a body segment with all the bones, joints and simulate muscles with varying levels of detail. With the ability to import the motion capture data and to perform forward and inverse dynamic simulations, it can calculate the force, displacement and torques on each of the body segments. This provides an important tool to calculate the forces acting on the L4 vertebrae without having to insert a pressure sensor into the subject. Because of this ability, the remodeling techniques can eventually be put to practical use.

The full body model generated by LifeMOD, along with the anthropometric data, treats the central torso containing the lumbar vertebrae as one segment (Figure 11). In order to calculate the loads on L4 vertebrae this model needed to be refined by adding individual vertebral bodies and intervertebral discs in the lumbar section. The components involved in creating a detailed lumbar section is shown in Figure 12. The adjacent segments, vertebral bodies, are connected to each other via spinal or the intervertebral discs. The intervertebral discs are kinematic constraint connecting adjacent vertebral bodies. These are implemented as a 6 degree of freedom force between any two given segments which include all the connective tissue between them. The disc is formed as a spring force. The values for stiffness and damping coefficient are used from the anthrompmetric data base in LifeMOD.

The spinal discs are connected to the vertebral segments at user specified landmark locations. A landmark is any significant location on a segments or bones in the model. It is an anatomical marker which represents position and orientation and moves with the segment to which it is attached to. The landmarks on the vertebral bodies were placed according to the locations specified in the literature [33]. The spinal discs connect to the inferior landmark on the vertebral body above them and superior landmark on the vertebral body below them (Figure 12). A detailed lumbar model was created in LifeMOD and incorporated in the full body model replacing the existing central torso section (Figure 13).

Figure 11.

Fullbody model with the central torso enlarged.

Figure 12.

Exploded view of the lumbar vertebrae with spinal disc

Figure 13.

Closeup of the detailed lumbar model created in LifeMOD.

2.4.3. Load calculation with MOCAP data

With the full body model refined to include the detailed lumbar section, Motion capture (MOCAP) data is imported in to LifeMOD. Along with the x,y and z coordinates of the marker for each time-step the MOCAP data file contains the subject specific data (age, height, weight and lengths of body segments). This information is used by LifeMOD to scale the musculoskeletal model and all its segments to match the subject data. Apart from the automatic scaling and positioning, manual tweaking to position the segments in the correct position by visual inspection of the landmark locations was necessary. The final body model along with the motion capture data is shown in Figure 14. The yellow spheres represent the marker agents imported through the motion capture data and the red spheres represent the locations on the body segments to which the marker agents are attached. Within LifeMOD the motion data (yellow spheres) and the marker locations on the body segment (red spheres) are linked via spring elements. The use of spring elements becomes clear in the next step. This helps to make adjustments for small differences in the human body model and the actual human subject.

Figure 14.

Body model along with the motion capture data.

With the motion capture data imported into LifeMOD, the next step is to perform an equilibrium analysis (or Settling run). In this process the model is fit to the data positions. The body segments are repositioned, while still maintaining the joint definitions, to match the position represented by the motion capture data. It is a dynamic analysis in which the data driven marker agents (yellow spheres) are fixed and the body segments are moved to find minimum energy configuration in the springs connecting the marker agents (red spheres). While recording the motion data, the markers placed on the body segment might not be exactly on the specified location. This means that, during equilibrium analysis the markers on the body segments will never overlap with the motion capture data points. The equilibrium analysis helps to reposition the model to align with the initial position of motion agent data.

Having the motion agents and the body segment markers agents connect via spring elements allows for minor differences between the locations of these two marker sets. The final body position after equilibrium analysis is shown in Figure 15.

Figure 15.

Body model after equilibrium analysis.

The next step is to perform the inverse dynamic simulation. During the inverse dynamic simulation, the joints and muscles are set as recording elements and the human body model is manipulated by the motion agents. The recording elements collect the angular motion and muscle lengthening and shortening data. After the inverse dynamic simulation, the motion agents are disabled and the recorded angular movements and the lengthening and shortening patterns are used in torque and force functions in joints and muscles. The joints and muscles are set to active from recording state. These active elements then drive the model along the curves recorded by the inverse-dynamics simulation as the PD-servo drivers replicate the joint and muscle data as closely as possible. At the end of the forward dynamic simulation the forces of the discs in all three directions (x, y and z axis) are calculated by LifeMOD.

The computed loads acting vertically down on L4 vertebra for two different activities (walking and rope jumping exercise) are shown in Figures 16 and 19. These values were comparable with load values for different activities in the literature [31, 32, 34]. The load values calculated from LifeMOD are used in the remodeling algorithm explained in the next section (2.5).

2.5. Adaptive bone remodeling

Many theoretical and computational models have been proposed to investigate and simulate the dynamic behavior of the bone [9, 35]. With high resolution imaging methodologies such as μCT becoming more accessible, the study of cancellous remodeling began to take in to account the influence of cellular activity in 2 and 3 dimensions [9, 25]. Simulation involving effect of metabolic bone formation deficit and the micro structural bone formation deficit were also developed [9, 36]. Effect of age-related changes on cancellous plate thickness and the number of plates in the cancellous bone showed that the number of plates per unit volume is a more significant predictor of bone loss than a decrease in the plate thickness [46].

2.5.1. Remodeling algorithm overview

Different mathematical control models of mechanical bone mass regulation have been proposed [37, 38]. These models assumed a continuous feedback loop between the maintenance of bone mass and local strain values in the tissue. The idea of functional adaptation of the cancellous structure regulated by the local activity of cells was proposed by Roux [39]. The local activity of the cells is governed by the mechanical stimuli. Mullender et al. [37] proposed that the osteocytes act as sensors of mechanical signals or mechanoreceptors because of their favorable architecture and distribution. The mathematical model uses the amplitude of strain energy density as the mechanical signal that the osteocyte senses. The osteocytes in turn send this signal to the bone surface. The signal intensity decreases exponentially with distance (Equation 2).

$$f_i(x)=e^{-(d_i(x)/D)}$$E2

where,

d_i= distance between osteocyte i and surface location x

D = Distance form osteocyte i where its effect

has reduced to e⁻¹ or 36.8%

The effective stimulus P(x, t) at any given surface location xis the sum of all the osteocyte signals within the regions of influence of nearby osteocytes (Equation 3) ([7, 37]).

$P(x,t)={\displaystyle \sum _{i=1}^N{f}_i(x){\mu }_i{R}_i(t)}$E3

where,

μ_i= Mechanosensitivity of the osteocyte iR_i(t) = Strain energy density rate

at osteocyte location i

The local change in the relative density mon the surface of cancellous bone tissue is expressed as [7],

$\frac{dm}{dt}=\{\begin{array}{c}\tau \{P(x,t)-k_{th}\}-r_{oc}ifP(x,t)>k_{th}\\ 0ifP(x,t)\le k_{th}\end{array}$E4

where k_this the threshold of bone formation, r_ocis the relative amount of bone resorbed by osteoblasts and τis the time constant regulating the rate of the process. This was chosen to be 1 making the simulation process measured in simulation time. The relative density at any location within the bone varies between 0 (for no bone) to 1 (fully mineralized bone). When the relative density value reaches zero, the osteocyte at that location disappears.

The resorption is regulated by the presence of micro cracks or by disuse [7, 8]. Resorption due to the presence of micro cracks is assumed to be constant,

$p(x,t)=d(spatiallyrandom)$E5

Resorption due to disuse is dependent on the strain and is given by the following equation, where a= 1.6 and c= 12.5 are constants.

$p(x,t)=\{\begin{array}{cc}c[a-P(x,t)]& ifP<a,\\ 0& ifP\ge a\end{array}$E6

Osteocytes respond to dynamic loads [40, 41]. The response of the osteocytes is measured in strain energy density (SED) rate. This SED rate can be determined by static finite element analysis using an equivalent static load of the dynamic loading signal. Consider a periodic loading signal σ_e(t), oscillating between 0 and a maximum amplitude of σ, given by Equation 7, which has a frequency of fHz (ω= 2π f).

$${\sigma }_e(t)=\frac{1}{2}\sigma (1+\cos \omega t)$$E7

Then the final form of equivalent static load (ESL) is given by Equation 8,

${\sigma }^{\prime }=2\sigma \sqrt{f}\Rightarrow ESL$E8

The previous studies have provided proof of the regulatory mechanism of the Wolff’s law. The idea proposed in this research is the use of the realistic mechanical loading data on the L4 vertebra (Section 2.4) to run the remodeling algorithm. The load data is then used in the remodeling algorithm to predict the modified cancellous structure. The schematic of the proposed method is shown in the Figure 16. The advantage of this is that the effect of different activities on the strength and structure of cancellous bone can be studied. Another advantage is by using the motion data of different everyday activities, a daily effective load value can be calculated and this in turn is used to drive the remodeling process to simulate the changes in the cancellous bone structure.

Figure 16.

Load acting in superior-inferior direction (vertically down) on L4 vertebra during walking.

2.5.2. Equivalent static load calculation (ESL)

The effective static load calculation of two different activities, walking and rope jumping exercise, are compared here. The load calculated from LifeMOD for walking motion of the subject is shown in Figure 16. The amplitude of the dynamic load oscillates between a minimum of 242.79N and a maximum of 866.04N, with a peak to peak amplitude of 623.25N. The first few data points were discarded as the model goes through a settling run in which the motion agent markers are repositioned to align with the motion agents. The last few data points were discarded to make sure none of the reflective markers on the subject were missing as a result of being close to or going out of the capture volume. The data retained is between the two vertical red lines at time 0.8 seconds and time 3.98 seconds.

The frequency of the load cycle is equal to the frequency of the subject walking at a self selected pace. Within the Vicon Nexus software two events of gait cycle, toe off and heel strike, are detected. A heel strike is defined as the instance the heel of the foot touches the ground and the toe off is defined as the instance when only the big toe of the reference foot is in contact with the ground. These two events, for both the feet, are shown in Figure 17. The toe offs are indicated by the red upward pointing arrows and the heel strikes are indicated by the black diamonds. Knowing the total duration of the capture data and the number of toe offs, the frequency of walking is determined. In this case, the first left foot toe off occurs at frame 331 and the last toe off occurs at frame 651 with the total number of frames in between being 338. With a total toe offs of 3 and at a capture rate of 100 frames/second, the time between toe-offs is 1.1267 seconds and the frequency equals to 0.8875Hz. With this information the effective static load can be calculated to be 1174.29N. This calculated load is used on the segment of cancellous bone in the finite element algorithm.

Figure 17.

Walking gait cycle showing heel strike and toe off.

The second activity was rope jumping exercise. The calculated load signal on L4 is shown in Figure 19. Once again the first and last few data points are discarded and signal from 1.732 seconds to 7.721 seconds is retained. The signal oscillates between a maximum amplitude of 1247N and a minimum amplitude of 145.9N. This has a peak to peak amplitude of 1101.1N for the dynamic signal. The frequency is determined similar to the method in walking gait. Between frames 110 and 804 for the first and last toe-off, there are 14 toe-offs (Figure 17). This gives 49.57 frames or 0.4957 seconds between toe-offs with a 100 frames per second capture rate. Thus the frequency of the load signal is 2.0173Hz. With this, the ESL is calculated to be 3127.82N. Once again, this calculated load is used on the segment of cancellous bone in finite element algorithm to calculate the stresses and strains within the bone segment.

Figure 18.

Rope jumping exercise showing heel strike and toe off.

Figure 19.

Load acting superior-inferior direction (vertically down) on L4 vertebra during rope jumping exercise.

2.5.3. Remodeling results

To calculate the stress and strain values within the bone material, a finite element solver specifically designed for biomechanical applications is used [44]. This is a three step process, sample preparation, finite element analysis and post processing of the data. These three steps are explained below.

A small volume of cancellous bone is extracted fromthe stack of μCT images. A 3D geometry of the sample is constructed. At this stage a tetrahedral mesh is created using a 3D tetrahedral finite element mesh generation algorithm [45]. A custom MATLAB code was written to assign material properties to the tetrahedral mesh and add a top and bottom rigid body plates with compressive loads on them [11]. Vertebral cancellous bone has an elastic modulus of 5,000 MPa and Poison’s ratio of 0.3 [8, 37]. The top and bottom of the bone segment are attached to rigid bodies whose translations are constrained along X and Y axes and rotation is constrained about X, Y and Z axes. A positive force (compressive) is applied on the rigid bodies equivalent to the calculated ESL from the previous section. The solver is set to perform a quasi-static analysis ignoring the inertial effects. After the solver completes the analysis, the stress and strain values at each element location are known. The strain energy density value is calculated first at the element location and then interpolated to find the value at each osteocyte locations.

The calculated strain energy density values mapped on the bone sample are shown in Figure 20. With the SED values calculated, along with the other parameters required for the remodeling algorithm explained in section 2.5.1, the bone goes through the process of increasing/decreasing the bone density values until the density value either goes above 1 or below 0.01. At that point, the voxel whose density goes below 0.01 is removed and new 6-neighbor voxels are added around the voxel whose density value goes above 1. Each time there is a change in the structure, new stress, strain and SED values are calculated.

Figure 20.

Calculated SED values mapped on the initial geometry of the bone sample.

This process is repeated until there is no further significant change in the structure between iterations. This is decided by comparing the quantifying parameters calculated at every 5th iteration (Figure 22, 24). The final structure of a bone sample for the corresponding action of walking is shown in Figure 21 overlapped with the SED values. Figure 23 shows the final geometry of the same bone sample for rope jumping exercise. The results show that the bone microarchitecture responds differently to different loads. In this initial test, only the vertical component (Z-axis) of the load is applied. Hence, the structure seems to evolve into vertical struts. This can be further refined by adding the horizontal components (X and Y-axes) of the load. In which case, the horizontal struts might be retained and/or modified to adapt to the loads in the horizontal directions.

Figure 21.

Final geometry of the bone structure (Applied load corresponds to walking).

Figure 22.

Variation in the quantifying parameters during remodeling as a result of walking exercise.

Figure 23.

Final geometry of the bone structure (Applied load corresponds to rope jumping exercise).

Figure 24.

Variation in the quantifying parameters during remodeling as a result of rope jumping exercise.

3. Conclusion and future work

Successful integration of the proposed method is demonstrated in this study. Effects of different physical activities on the adaptation of the cancellous bone structure and their quantifying parameters are shown. The results suggest that the different loadings corresponding to different activities have varying effect on the microarchitecture of the bone. Thus, the loading conditions on the bone can be changed to produce desired results, i.e., stronger or at least not a weak structure. In general, this method has applications in physiotherapy, assessing/avoiding injuries in various sports, and studying the effects of different types of bone implants over a period of time.

Although the proposed remodeling cycle is successfully integrated and implemented in this study, there are still rooms for improvement. Some of the limitations of this study which would be the focus of the future work are, 1) the μCT scans of the vertebra are not in vivo, 2) effect of only one motion is included in the simulation which can be refined by including the daily loading activity, 3) Loading on the only the vertical direction is considered which may be extended to include all the three forces. Some of the other possibilities for improving the performance would be to look into the use of Graphical Processing Units (GPU) for computation. Use of GPU may be more cost efficient in a clinical environment. Based on the results, this study suggests the above mentioned to be addressed in the future work.

Acknowledgement

The authors would like to thank Dr. Kathleen Issen and Dr. Laurel Kuxhaus for their valuable advice and constructional critics and Nicole Corbiere for helping with μCT scanning validation.