Calculated parameters of the metal foam sample.
Vertebral compression fractures are very common low-bone-density related conditions among the elderly and have a significant adverse health effect . 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 . The fractures affects up to 25% of postmenopausal women . Men are equally affected by this, but because of greater peak bone mass, it occurs about 10 years later than women . 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 . 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 . 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
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 . 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
An important step in the analysis of
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).
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
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.  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,
Find the bounday voxel
Find the critical points in the set of boundary voxels
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
Remove the non-critical points from the set of bondary points
Repeate steps 1-3 untill there are no more boundary points to remove
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 .
Next, the algorithm was tested on a CAD model with plates in it and
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.
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 . 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 .
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 . A prominent change of plate-like structures to rod-like structures was shown by Ding et al.  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. . 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 . Recent advances in
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.  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 . Tetrahedral mesh generation algorithm  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 . According to the volume-based approach, the local thickness (
Figure 8 illustrates this concept. The blue line is the object boundary, red line is the mid-axis of the object. If
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
The quantifying parameters were measured on this CAD after rotating it around the
The third structure used for verification was a metal foam sample scanned at a resolution of 28.7
|Parameters||Calculated values||Measured values||% Error|
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 . 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 .
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 . 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.
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 . 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).
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.
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.
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
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 . The local activity of the cells is governed by the mechanical stimuli. Mullender et al.  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).
The local change in the relative density
Resorption due to disuse is dependent on the strain and is given by the following equation, where
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
Then the final form of equivalent static load (ESL) is given by Equation 8,
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.
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.
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.
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 . 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
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.
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.
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
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.