The extended SIBA algorithm is based on the remodeling sequence first described by Frost [7], which consists of four phases: resorption; reversal; formation; and quiescence. Thomsen et al. [1] used an iterative approach for their stochastic model and the time interval between each iteration was one day. In our approach [8], one iteration would include all four steps, where bone resorption (osteoclastic activity) and bone formation (osteoblastic activity) would be modeled for each iteration. Nevertheless, the timeframe for the total remodeling cycle was chosen to equal 197 days to meet the average remodeling cycle time from the Thomsen study. Furthermore, the algorithm is based on the assumption that for normal bone osteoclastic and osteoblastic activity are in balance resulting in an almost stable bone mass and bone structure. After the onset of menopause, which was assumed to be at the age of 50 years, the balance would become negative, i.e. more bone is resorbed than formed resulting in a net loss of bone mass. No further information was needed to run the models. The changes in architecture, we observed, were therefore only the result of the loss in bone mass over the years, which means that the architecture as it was present before menopause already predetermined the outcome in the simulated bone atrophy with a given, constant loss of bone. In the simulation, osteoclastic activity is modeled as an exponential function, i.e. the probability for bone being resorbed at a certain location in the bone matrix has a Gaussian profile. In practical terms this means that the algorithm, applied to the digitized trabecular bone structure, is based on constrained Gauss filtration of binary data sets, i.e. only bone and nonbone will be differentiated. The effect of the Gaussian filter is mainly controlled by the filter width and the constrained filter support. Constrained, in this sense, means that the discrete convolution of theGaussian and the binary data set is computed only for a limited, finite filter support. In our case, we used a filter support of 2, which means that a newly computed value only depended on the sum of the weighted values of the point's direct neighborhood (5 x 5 x 5). This also means that an osteoclast can penetrate into the bone only two voxels per cycle which corresponds to a perforation depth ranging between 28 and 48 µm depending on the neighborhood configuration. Since the original structures were based on binary data, every data point in the volume has either a value of one or zero prior to the filtration. The constrained Gauss filtration will result in a grayscale data set with values between zero and one. The more data points are 'set' in the weighted neighborhood the higher the resulting grayscale values and vice versa. Therefore, the application of the Gaussian filter introduces a small gradient on the surface of the structure, whereas values in the center of the structure, totally surrounded by mineralized bone, will be unchanged. The total amount of bone mass in a local neighborhood (sum of weighted voxels) is the same before and after the filtration but the threedimensional distribution of the mass has changed according to the local topology. In our simulation, the four phases in the remodeling cycle are fully synchronized for all cell activities in the structure, which means that in a first step bone is resorbed only and then in a second step bone is formed only. There is no temporal overlay of the osteoclastic and osteoblastic activity, which in reality is not the case and has therefore to be taken as a current limitation of the algorithm. Additionally, the amount of surface area involved in a remodeling cycle was not controlled, so that in the simulation the whole surface was affected by wither osteoclastic or osteoblastic cell activity, which again in reality is most probably not the case. The amount of bone that is removed or added in one iteration step is controlled by setting a global threshold for the filtered data, which must be between zero and one. In our simulation we found, that a threshold of 0.5 (all voxels with a grayscale value higher or equal to 0.5 are included) would allow to maintain the total amount of mass. This was defined as the 100% efficiency level for the osteoblastic activity which means that the osteoblast were able to form as much bone as has been resorbed by the osteoclasts. A threshold of 1.0 on the other hand resulted in the removal of all elements with at least one zero voxel in their local 5 x 5 x 5 neighborhood, which would correspond to a 0% efficiency level. The relation between the threshold and the efficiency level can be expressed as a second order polynomial, which means that initially the efficiency is high and is dropping down quickly with increasing threshold.
Whereas the original implementation distinguished only between pre- and post-menopause [8], where changes were assumed to happen only after the onset of menopause, the current implementation includes three distinctively different stages, namely the pre-, peri- and post-menopause. Pre-menopause is the period of time after bone growth ceases and before menopause begins. Although bone loss in this stage is minimal, it is important to model the structural changes prior menopause. Without its inclusion, the changes that occur will be underestimated. The peri-menopausal stage is a period of about five years after the onset of menopause where the bone may actually undergo the most changes [9,10]. Experimental values from the literature were used to determine the values for osteoclastic perforation depth, the osteoblastic efficiency level and activation frequency. Activation frequency (AF) is defined as the rate at which bone multicellular units (BMU) are formed. Together with the individual cellular function rates at the BMU, it determines the percentage of bone turnover [11]. AF was not part of the original algorithm [8], where the length of one iteration was defined by the length of the remodeling cycle (197 days). In the present simulation, the length was simply defined as 1/AF under the assumption that each BMU is activated in that time. Normal activation frequency is 1/1096 days [1], which corresponds to an iteration length of approximately 3 years. In peri-menopause activation frequency is doubled [12], which results in a cycle time of 1.5 years.
In order to establish a more realistic time frame and an estimation for the values of simulation input parameters, the experimental postmenopausal group was used as a control group. The parameters used for the three stages may differ for each individual but were chosen to be identical for all specimens within a specific stage with an osteoclastic penetration depth of approximately 30 µm (sigma=1.1, support=2), and an osteoclastic efficiency level of 99.7% (resulting in a decrease in volume density of 0.34% per year [13]) and an AF of 1/1096 days for the pre-menopausal phase. The number of interations needed to reach the age of 50 years, the assumed onset of menopause, varied between zero and nine depending on the individual age of the donor bone. The average bone loss that occured was -5.6% over an average of 13.4 years of pre-menopause, which was consistent with the values found in the literature for that time period [13]. In peri-menopause the osteoclastic penetration depth was increased by 40% changing the sigma from 1.1 to 1.5, the osteoclastic efficiency level was reduced to 95% and AF was doubled to 1/1096 days resulting in a total of four iterations needed to simulate six years of peri-menopause. Finally, in post-menopause, osteoclastic penetration depth and AF were set to the values for the pre-menopausal phase with a sigma of 1.1 and an AF of 1/1096 days respectively. These values produced a mean bone volume density at the 77-year-point of the simulation almost identical to the experimental post-menopausal group. Post-menopause was simulated with a total of eight iterations concluding at the age of 80 years.
It is clear that this approach is only one of many scenarios of what could have happened to these premenopausal bone specimens in time. Important for us was to see whether the outcome of the simulation would result in reasonable and realistic structures if compared to the postmenopausal group. The validation of the algorithm is based on both a qualitative (visualization) and quantitative (morphometry) comparison between those two groups. For the threedimensional visualization, an extended marching cubes algorithm was used [2]. The algorithm, which is based on a divideandconquer approach, triangulates the surface of any given voxel array in a fast and secure way. It does not only allow perfect triangulation but also smooths the surface effectively.
References: