CC BY
17
MSC 62P10
DOI: 10.14529/mmpl80204
A PARAMETRIC STOCHASTIC MODEL OF BONE GEOMETRY
V.I. Zalyapin1, Yu.S. Timofeev2, E.A. Shishkina2
^outh Ural State University, Chelyabinsk, Russian Federation, 2
E-mail: zaliapinviiQsusu.ru, swimmer_86iQmail.ru, lenaiQurcrm.ru
The aim of the present study is to develop a parametric bone modelling algorithm which takes into account bone microarchitecture. This approach allows to generate hematopoietic bone segment phantoms based on literature-derived micro- and macro dimensions. We propose a method for subdividing bones into small segments which can be described by simple geometric shapes filled with a stochastically generated rod-like model of the trabecular structure with appropriate voxel resolution. This approach avoids the disadvantages of non-parametric individual modelling based on computer tomography scans. The parametric method allows the simulation of individual variability in bone-specific dimensions. The model presented in this paper will be used to describe the geometry of hematopoietic sites, which in turn will serve as a basis for calculating the doses of irradiation of the hematopoietic cells of the bone marrow from the incorporated beta-emitters.
Keywords: micro- and macro- structure of the trabecular bone; stochastic modelling; voxeliza,tion.
Introduction
The bones of the human skeleton are composed of spongiosa covered with a compact bone layer, or cortical bone (Fig. 1). Spongiosa is a porous structure with friable lying bone trabeculae permeating bone marrow the soft tissue of the internal bone cavity, which includes hematopoietic cells (active bone marrow).
Fig. 1. Outer cortex and internal spongy structure (trabecular bone) of the human lumbar vertebra
Skeleton geometry modelling is widely used in biomedical research, for example, in studying the mechanical properties of bone in traumatology and implantation [1]; in studies of relationships between bone morphology and different diseases (osteoarthritis, osteoporosis) [2]; in radiation dosimetry [3 5], and so on. We are interested in the latter problem.
Different tasks require different levels of geometric details and different degrees of realism. For example, very realistic individual models of a patients bone shape are created for medical examination with gamma irradiation (computer tomography, cyber knife.
SPECT systems). Hovever, there is no need to build a material microstructure model; both cortical and sponge bones can be considered as homogeneous structures in terms of photon transport [3]. Bone dosimetry for alpha emitters doesn't require a description of the bone shape (due to short pathlength of a-particles). The dimensions of a bone segment model are usually not larger than 2 mm [4]. However, it is extremely important to describe the trabecular structure at the sites of radionuclide localization (for instance, trabecula surfaces and adjacent intertrabecular space) as accurately as possible.
Dose estimation in active marrow exposed to bone-seeking beta-emitters, such as 89Sr or 90Sr/90Y (energies of electron emission are 0 - 1,5 MeV and 0 - 2,4 MeV, respectively), is an important task in bone dosimetry. Monte Carlo simulations of electron - photon transport to calculate the active marrow doses are based on the geometrical modelling of bone structures [4,5,19]. The description of spongy bone microarchitecture in the model is important due to high probability of low-energy electron emission. The mean free paths of
90 90
in the range of 1-3 mm, depending on the density of the bone media. These dimensions are comparable to the thickness of the cortical layer of skeletal bones. Therefore, the evaluation of the contribution of contaminated cortical bone into bone marrow dose relies on the accurate description of cortical thickness. The maximum free path (for electrons yielding with maximum energy) in a spongy bone can reach 5-9 mm. These dimensions are comparable to the linear dimensions of many bones of the skeleton. Thus, the model geometry should consist of accurate descriptions of the fine structure of spongiosa and cortical bone thickness, as well as descriptions of macro-dimensions of the bone.
New computer tomography (CT)-based methods are widely applied to developed computational phantoms [5], among which there are phantoms designed for bone dosimetry of beta-emitting radionuclides [6-8]. Such phantoms were developed using a combination of micro- and macro- CT images of human bones. The advantage of this CT-based method lies in the high realism of the complex bone shape description, as well as the possibility of an adequate description of bone-specific microstructure. However, the method has a number of disadvantages, such as:
• The computations are laborious and expensive.
•
exposure. Cadavers are commonly used, associated with organizational difficulties.
estimation of uncertainties associated with the variability of human anatomy.
•
of the imitated structural components to allow accurate simulation of electron-photon transport. In practice, the voxel resolution of the model is limited by computational power. For example, the voxel sizes of the models describing large bone segments [6-8] were restricted to 100 - 200 fim, despite the fact that the thickness of the trabeculae (the modeled structural elements) ranges from 40 to 400 fim. Inadequate voxelization can lead to errors in dose calculations. In theory, it is possible to change the voxelization when increased computing power is available, but in practice this is difficult to implement, since it would be necessary to repeat the analysis and alignment of CT and ^CT images.
1mm) is insufficient to determine cortical thickness for most of bones (for example, about 0,4mm in vertebral bodies [9].
Moreover, high individual variability of bone shapes and macro-dimensions negates the advantages of the high realism of the model of a single (or a population-average) body.
The aim of the present study is to develop an algorithm of parametric bone modelling, which allows the generation of phantoms of hematopoietic bone segments based on literature-derived micro- and macro dimensions.
We propose an approach that permits easy subdivision of bones into small segments which can be described by simple-shape geometric figures with appropriate voxel resolution. The parameters of the proposed cadaver-free model can be obtained from publications on morphometry and hystomorphometry widely presented in the biomedical literature. The proposed cadaver-free method avoids the disadvantages of non-parametric individual modelling based on CT scans. The parametric approach allows the simulation of individual variability of bone-specific dimensions.
1. Basic Assumptions
We will proceed from the following assumptions.
1. According to preliminary estimates [10]. variations in the geometric shape of bones with similar linear dimensions do not introduce a significant error (<5%) in the calculation of doses. Therefore, it is possible to describe the shape of the bone with simple geometric figures (stylized phantoms).
2. The cortical layer is considered as a homogeneous and isotropic substance located between two surfaces, one of which is the outer boundary of the stylized phantom, and the other is separated from it inside the phantom by a distance equal to the cortical thickness. Cortical thickness (Cb.Th) is a parameter of the model, which can be estimated from literature data.
3. The trabecular bone can be described by an isotropic stochastically deformed network structure, with beams randomly varying in width. This model can parametrized using standard literature-derived indicators of bone microarchitecture, such as:
• mean trabecular thickness (Tb. Th) and its intra-specimen variability;
•
intra-specimen variability; and
and the range of its admissible values.
4. The generated analytical model must be voxelized at a given resolution.
2. The Bone Shape Model
The following geometric objects are used as the stylized phantoms designed to describe the shape of bones:
Fig. 2. Stylized phantoms
The parameters of these figures are known from the literature and are specified by the
For example, a fragment of a rib body can be described by a rectangular parallelepiped, as shown in Fig. 3. Standardized morphometric indices, such as h and w (see. for example. [11]) are averaged over the entire length of the edge and are attributed to the height and width of the model parallelepiped. The length of the segment of the edge is chosen on the order of 3 lengths of the maximum range of electrons (in the example presented. 3 cm). The average value of the cortical thickness is also taken on the basis of the literature data [12]. Similarly, the neck of the hip is modelled by a cylinder, the head of the hip and humerus is half-ellipsoidal, the wings of the sacrum are triangular, the region of the spines of the femur and humerus is truncated cones, the acetabulum is a double cylinder, etc.
Fig. 3. An example of creating a stylized rib phantom
The main microparameter at this stage is Cb.Th. The stylized model of the shape of the bone and its cortical layer encompasses a region of space filled with a trabecular structure. The trabecular structure of the bone is simulated by immersing the phantom in an enclosing rectangular parallelepiped (carcass), and drawing the main parameters, as discussed below.
3. A Stochastic Model of the Trabecular Bone Structure
The microstructure of the trabecular bone is modelled in four stages: construction of the framework, randomization of the framework nodes, stochastic modelling of rod-like bone structures, and calibration of the model.
3.1. Construction of the Frame
The main microparameters at this stage are trabecular separation (Tb.Sp) and trabecular thickness (Tb.Th). In a simple initial frame with cubic cells Tb.Sp will be equal to the length of the edges of cells, or conditional trabeculae. According to data from the literature [13 15], Tb.Sp can be considered a normally distributed value, while Tb.Th has a lognormal distribution.
We take the correct frame in R3 composed of rectangular parallelepipeds with edges parallel to the coordinate axes K = {(ti,t2,t3) : ai < ti < bi, i = 1, 2, 3.} to construct the framework (Fig. 4. left). The dimensions of the frame are chosen so that the figures of the stylized phantoms (F) are completely immersed in them: VK>VF.
Fig. 4. The framework for microstructure modelling (left) and grid cell notation (right)
Dividing abscissa (¿i), ordinate (t2) and applicate (¿3) into m, n and l parts, respectively, the frame nodes can be denoted as follows:
a1
<ti < ... < tn
bi, a2
¿2 < ¿2 < ■■■ < ¿2
bo
a3
¿2 < ¿3 < ■■■ < ¿3
b3
where each of the frame nodes is denoted by the ordered number triplet (i,j,k), corresponding to the node coordinate Mijk = (¿1, ). The notation Ml]k = (¿{,¿2,^)) will always refer to a grid cell as shown in Fig. 4 (right), so that grid cell numbers range (n, m, l)
The parameters of the framework partitioning, (i.e., the values A^ = ¿S+1 — ¿s, i = 1, 2, 3, s = 0,1, 2,Ni, N{ = n, N2 = m, N3 = l,) are determined by the characteristic distances between the trabeculae in the horizontal, vertical and frontal directions.
As a first approximation, we can assume that the values AtS = ¿S+1 — ¿S are sufficiently close to Tb.Sp the characteristic distance between the trabeculae, so that the horizontal, vertical and frontal dimensions of the frame parallelepiped are multiples of this quantity.
3.2. Node Randomization
/6 1
Let £ = < £2denote a Gaussian random variable with the zero vector of averages
I £3
and the covariance matrix K%
(al ki2 ik 3
K = I ki2 to to k23
3 k23 q2
It should be noted that we assume the node permutations are mutually independent events. Each permuted node Mijk can be expressed in terms of its initial state M = Mijk according
to the following rule (Fig. 5): M = M + £
M
tj = j + £Sjk, s = 1, 2, 3
Fig. 5. Transformation of the framework elements and the volume calculation
Thus, the transformed framework is based on the natural transformation of the original edges in accordance with the transformed node positions. For example, an edge initially linking the nodes Mijk and Mi+ijk will turn into an edge linking the mapped nodes {MIijk and Mi+ijk). The randomization parameters (that is, the values a2, kij) are determined by the scattering and correlation relationships of the value Tb.Sp1. As the first approximation, we can assume that a2 = a2, kij = 0.
3.3. Parameters of the Modified Framework
If A(t^, ttA1) and B{tB, tB, tB) are arbitrary adjacent nodes, then the length of corresponding edge can be expressed as follows:
\AB\ = ,j(tf + tf - tA - tf)2 + (tB + {B - tA - &)2 + (tB + - tA -
For example, the length of transformed abscissa-directed (frontal) edges (MijkMi+1jk ^ MjkMi+ijk) can be calculated as
\Ml3k M%
i+ljk
\ = 1+1+ei+1 ,jk - ej -1\)2 + (c2+1 j k - ej )2 + (e^1 ,jk - j )2 •
The length of transformed ordinate-directed (lateral) edges, MijkM^+ik M^kMj+ik, and applicate-directed (vertical) edges, MijkMijk — MjkMlijk+i edges can be calculated similarly.
Let Atpss ,ps = 1, 2,...,Ns,s = 1, 2, 3 denote the distance between the parallel edges of the parent framework, where s = 1 corresponds to the frontal direction; s = 2 to the
s=3
example, in the frontal direction),
A£
A£i
Л£2 Л£з
¿i+1jk ¿ijk e1 - e1 ¿i+1jk i-ijk e2 - e2 ei+1jk _ eijk S3 S3
is the Gaussian [13] random variable with the zero vector of averages and the covariance matrix
K,
A?
(
\
22 a1(i+1jk) + °1(ijk)
k12
k13
k12 2
a2(i+1jk) k23
+ a,
2(ijk)
kA? kA?
k23
\
a.
2
3(i+1,jk)
+ a
2
3(ijk)
where kpq
ki+1jk I kijk kpd + kpd ,
p,q = 1, 2, 3. A similar expression can be easily written for other two directions. If the permutation scattering is significantly smaller than the typical
1Here it is useful to bear in mind the following circumstance: the average value of the intertrabecular space is very unstable ([11]) and can vary within fairly wide limits. Therefore, the adjustment and selection of the values for the parameter Tb.Sp take into account the behaviour of the ratio of the volume of bone tissue (BV) to the total bone volume (TV). The desired value of this ratio is set by the user as a range of possible values before starting the simulation. If the resulting value of the BV / TV ratio differs from the desired value slightly - the parameters Tb.Sp are accepted. Otherwise, they are adjusted so that the BV / TV ratio falls within the desired range.
2
size of parent framework structures, << (Atj )2, then a first approximation of the
average distance between edges of the transformed framework can be calculated as2 :
- ^ >=*+(^ )2+(j )2+(c )2+(j )2+()2+(j )2 •
This relationship allows us to specify the parameters of the original framework in order to obtain the desired parameters of the modified one.
3.4. Cell Volume of the Transformed Frame
Let A1A2A3A4B1B2B3B4 be the nodes of the transformed frame for a given cell (ijk) (Fig. 5)3. The cell can be subdivided into two prisms, A1A2A4B1B2B4 and A2A3A4B2B3B4, and each of the prisms can be further subdivided into 3 triangular pyramids as shown in Fig.5 to simplify the volume calculation.
Pyramid volume can be calculated as V = 11(a, b, c)\, where ( a, b, c) is the triple scalar product of the reference vectors of a pyramid.
Assuming (as above) &2(ijk) << (Atj)2 and considering only those elements of triple scalar product whose degree of smallness does not exceeded vjj)' ^he v°lume °f the pyramid can be calculated as follows4:
v = 1 ((ti - ti~i)(tj2 - t2~l)(tk - tk~i) + (tr1"-1 - e-lj-lk~l№ - t^m -111)+
+(j-1 - ei-1k-1№ - ti-1)(tk - tk-1) + j - %-1k-1)(A - ti-i)(tj - tj-1)).
3.5. Modelling of Stem-Like Bone Structures
We assume that a single trabecula is located on every edge of the modified frame (Fig. 6). The distribution of trabecular thicknesses within most bone samples is somewhat asymmetric and can be approximated by a lognormal distribution [14,15]. In each modified node Mjk we generate a random variable Rijk with given characteristics M(Rijk) = R, a2(Ri,jk) = aj. If Mi+1 jk is a node adjacent to the node Mlijk, then the trabecula will be modeled as a truncated cone with the axis MjkMi+1jk and spherical surfaces of radii Rjk and Ri+1jk respectively at the bases.
The trabecula volume (Fig. 6) can be calculated in 4 steps:
R1 , R2
the length is L)\
V(L,R1 ,Rj) = 3nL(Rj + RR + Rj), L = ^lj - (R+j - Rijk)j - (R1 - Rj)j,
D Ri+1jk\/12 - (Ri+1jk - Rijk)2 D Rijk\J12 - (Ri+1jk - Rijk)2
R1 = -j-, R2
l ' 2 l
2 Hereinafter, all equations are presented for frontal-directed perturbations. Lateral- and vertical-directed perturbations can be described by direct analogy.
3The node B2 in the current notation corresponds to Mjk
4It should be noted that the average volume of randomly transformed frameworks is practically equal to the volume of the parent framework. In contrast, the average edge lengths become larger after transformation.
Fig. 6. Trabecula volume calculation
2. Calculation of the volume of the spherical segments: VSi = nH2(R — |H) where H can be computed for two nodes using the parameters shown in the Fig.6. For the left node (Rijk < Ri+ijk) H2 = Rijk + u2, u2 = \JRjk — R, and f°r the right node
{Rijk < Ri+jk) :Hi = Ri+ijk — uij ui = \JR2+ijk — r2•
3. Subtraction of the volume of the supplementing spherical segments V (L,Ri, R2)—VSi.
4. Calculation of the overall trabecular volume as a sum of volumes of conoidal structures (after subtraction of the supplementing spherical segments) and spherical structures at the nodes of the transformed grid.
4. Model Calibration
It should be noted that Tb.Sp is an uncertain parameter and if adopted in the rodlike model it may lead to incorrect estimation of the spongioza density. In other words, this model can only be considered the first approximation. It is therefore important to fit a criterion of adequacy, where an acceptable range for the volume ratio BV/TV is adopted (see above footnote 1). in order to obtain model volume ratios corresponding to values observed in reality. For this purpose, we have developed a multi-step procedure for calibrating Tb.Sp values. The calibration is performed by changing the average distance between nodes at each step until the model ratio (BV/TV)i is equal to the calibration value (BV/TV)caiibr specified by the user before starting the simulation.
The iterative calibration algorithm is as follows. An additional coefficient (Tb.Sp.Factor) is introduced into the formula for calculating the distance between nodes:
Ati = rs+1 - rs = Tb.Sp.Factor • Tb.Sp. The perturbation vector is also multiplied by this coefficient:
I Afi+1 \ I AG \
A£+i = | A£2+1 I = Tb.Sp.F actor ■ | I .
V Af3+1 / V A ay
At the first calibration step, Tb.Sp.Factor is assumed to be equal to 1. We model the cellular structure, obtain the model ratio (BV/TV)j at the ¿-th step, and compare it to the calibration value (BV/TV)Calibr. If they do not coincide, the new value Tb.Sp.Factor-is calculated according to the following formula:
<Tb.sp№+i. (B/V^+1).
where i is the index of the calibration step.
When the DV/TV ratio calculated for the current model falls within the admissible range, the obtained model is considered calibrated, and the associated characteristics of the trabecular structure can be used for further calculations.
5. Model Voxelization and Cortical Layer Description
The final simulation result should be presented in the voxel form for further use in dosimetric calculations and visualization. The subdivision into elementary volumes (voxels) is carried out by a grid of planes located at equal distances from each other and perpendicular to each of the coordinate axes. The distance between the planes corresponds to voxel resolution.
We assume that each voxel contains one of the four homogeneous fillings:
(1) bone marrow;
(2) bone substance of trabeculae;
(3) bone substance of the cortical layer;
(4) empty space.
Thus, the simulated body is a Boolean four-dimensional array of voxels, represented by three coordinates corresponding to one of the corners of the cube-voxel, and a value from 1 to 4 corresponding to filling the voxel.
Unlike the analytically modeled trabecular structure and surfaces corresponding to the shape of the bone, the cortical layer model is created during the voxelization process. For a voxel partition in the volume of the space on which the stochastic model of the trabeculae (VK) is defined, we construct a parallelepiped (P), in which the modeled figure (F) is fully contained and for which the lengths of its edges are multiples of the voxel resolution VF < VP < VK. Here VF,VP,VK are the volumes of the respective figures. We determine whether a voxel belongs to a particular structure based on the central point of the voxel via the following steps:
1. If the coordinate of the center point of the voxel does not belong to to the space bounded by the surfaces of the figure F, then the voxel is assigned material 4.
2. Otherwise, we determine how far the center point is located from the nearest face of F. If the distance is less than or equal to the thickness of the cortical bone, then material 3 is assigned to the voxel.
3. If the previous steps do not yield a result, material 2 is selected if the central point belongs to the space corresponding to the elements of the rod-like structure of the trabeculae; and material 1 if it does not.
Fig. 7 shows the result of the voxelization of the rib model, previously described as an example (see Fig. 3). In this example, no empty space was generated. Bone structures (materials 3 and 2) are shown in black. Bone marrow (material 1) is shown in white.
Fig. 7. Visualization of the horizontal section of a rib segment
6. Simulation of Individual Variability
Simulation of individual variability is achieved by perturbing the input parameters of the model5. For this purpose, we construct a number of models with randomly generated main parameters based on the population-average macro-dimensions, Tb. Th, calibrated Tb.Sp and their person-to-person variability (namely, user-defined literature-based standard deviations). The calibration factor (Tb.Sp.Factor) for all of these models is the same. Each of the resulting models is subjected to verification; here the calculated BV / TV ratio must fall into the admissible range.
The drawing of micro- and macro-parameters is carried out with lognormal and normal approaches respectively, within 90% of the confidence interval. The values must agree as follows:
1- If are the random variables describing the macro parameters of the model, and ^ takes a value x1; then the value £2 = x2 agrees with the value ^ = x1 if F^2 (x2) = F^ (x1) Here F^ (x) is the cumulative distribution function of the random value The value can be easily obtained as x2 = F?"1(F?1 (x1)).
2- If are the random variables describing the micro parameters (Tb.Th and Tb.Sp) of the model, and ^ takes a value x1; then value £2 = x2 agrees with the value ^ = x1 if F^2 (x2) = 1 — F^ (x1) Here F^(x) is the cumulative distribution function of the random value The value can be easily obtained as x2 = F—1 (1 — F^ 1(x1)).
The Box - Muller [20] algorithm was used to draw the normal and lognormal random variables.
Conclusion
Originally, parametric stochastic model of bone geometry was developed to be applied for the radiation risk assessments in the South Urals. Radioactive contamination in the South Urals in the 1950s led to radiation exposure to people who were subsequently included in epidemiological cohorts to investigate radiation risk. Among the radionuclides which contaminated the environment, affecting the residents along the Techa riverside and in the East Urals Radioactive Trace, Strontium isotopes were the most dangerous for bone marrow exposure. For example, a statistically significant relationship between the radiation dose and leukemia cases was found in the Techa River cohort [17]. Thus, the assessment of the doses from the incorporated 89'90Sr, aimed at improving the Dosimetric System of the River Techa (TRDS) [18], is of paramount importance in the dosimetric studies currently under way in the JCCRER 1.1 project.
The model presented in this paper will be used to describe the geometry of hematopoietic sites, which in turn will serve as the basis for calculating the doses of irradiation of the hematopoietic cells of the bone marrow.
Currently, work is under way to implement the presented algorithm in the form of a computer programme with an interface designed for radiation biophysicists. The output will be compatible with the existing programme for modelling the transport of emissions.
5Note that when constructing a set of models, the macro parameters of the bone (linear dimensions) are subject to the normal distribution. Macro-parameters are assumed to be connected. (If one of the parameters turned out to be higher than the average, then for the others we expect values greater than the average.) Cortical thickness is an independent parameter. When modelling the variants of the spongy structure, its parameters (trabecular thickness and trabecular space) are also related: the greater the thickness of the trabeculae, the less the trabecular separation.
Besides incident-specific study, the model will allow estimating dosimetric uncertainties due to individual variability of bone geometry for the first time. Moreover, the results obtained with the model will be useful for verification of complex CT- and ^CT- based models of ICEP [6-8].
The work was funded by the Federal Medical-Biological Agency of the Russian Federation and the US Department of Energy's Office of International Health Programs in the framework of joint US-Russian Project 1.1.
References
1. Cox L.G.E., van Rietbergen B., van Donkelaar C.C., Ito K. Analysis of Bone Architecture Sensitivity for Changes in Mechanical Loading, Cellular Activity, Mechanotransduction and Tissue Properties. Biomechanics and Modeling in Mechanobiology, 2011, vol. 10, no. 2, pp. 701-712. DOI: 10.1007/sl0237-010-0267-x
2. Vaughan T.J., Voisin M., Niebur G.L., McNamara L.M. Multiscale Modeling of Trabecular Bone Marrow: Understanding the Micromechanical Environment of Mesenchymal Stem Cells During Osteoporosis. Journal of Biomechanical Engineering, vol. 137, no. 1, p. 011003. DOI: 10.1115/1.4028986
3. McParland B.J. Nuclear Medicine Radiation Dosimetry: Advanced Theoretical Principles. London, Springer Science and Business Media, 2010. DOI: 10.1007/978-1-84882-126-2
4. Dant J.T., Richardson R.B., Nie L.H. Monte Carlo Simulation of Age-Dependent Radiation Dose from Alpha- and Beta-Emitting Radionuclides to Critical Trabecular Bone and Bone Marrow Targets. Physics in Medicine and Biology, 2013, vol. 58, no. 10, pp. 3301-3319. DOI: 10.1088/0031-9155/58/10/3301
5. Bolch W.E., Lee C., Wayson M., Johnson P. Hybrid Computational Phantoms for Medical Dose Reconstruction. Radiation and Environmental Biophysics, 2010, vol. 49, no. 2, pp. 155-168. DOI: 10.1007/s00411-009-0260-x
6. Pafundi D., Rajon D., Jokisch D., Lee C., Bolch W. An Image-Based Skeletal Dosimetry Model for the ICRP Reference Newborn - Internal Electron Sources. Physics in Medicine and Biology, 2010, no. 7, vol. 55, pp. 1785-1814. DOI: 10.1088/0031-9155/55/7/002
7. Hough M., Johnson P., Rajon D., Jokisch D., Lee C., Bolch W. An Image-Based Skeletal Dosimetry Model for the ICRP Reference Adult Male - Internal Electron Sources. Physics in Medicine and Biology, 2011, vol. 56, no. 8, pp. 2309-2346. DOI: 10.1088/0031-9155/56/8/001
8. O'Reilly Sh.E., DeWeese L.S., Maynard M.R., Rajon D., Wayson M.B., Marshall E.L., Bolch W. An Image-Based Skeletal Dosimetry Model for the ICRP Reference Adult Female - Internal Electron Sources. Physics in Medicine and Biology, 2016, vol. 61, no. 24, pp. 8794-8824. DOI: 10.1088/1361-6560/61/24/8794
9. Fazzalari N.L., Parkinson I.H., Fogg Q.A., Sutton-Smith P. Antero-Postero Differences in Cortical Thickness and Cortical Porosity of T12 to L5 Vertebral Bodies. Joint Bone Spine, 2006, vol. 73, no. 3, pp. 293-297. DOI: 10.1016/j.jbspin.2005.03.023
10. Volchkova A.Yu.. Shishkina E.A., Tolstykh E.I. Effect of the Shape of the Bone on the Power of the Doses in the Bone Marrow of a Person. Materials of the VI International Scientific and Practical Conference, Chelyabinsk, November 8-9, 2016, pp. 36-41. (in Russian)
11. Roberts S., Chen P. On Some Geometric Properties of Human Ribs. Journal of Biomechanics, 1970, vol. 3, pp. 527-545. DOI: 10.1016/0021-9290(70)90037-0
12. Li Z., Kindig M.W., Subit D., Kent R.W. Influence of Mesh Density, Cortical Thickness and Material Properties on Human Rib Fracture Prediction. Medical Engineering and Physics, 2010, vol. 32, no. 9, pp. 998-1008. DOI: 10.1016/j.medengphy.2010.06.015
13. Sutton-Smith P., Beard H., Fazzalari N.L. Quantitative Backscattered Electron Microscopy of Normal and Osteoporotic Bone. Proceedings of the Annual Scientific Meeting of the Australian and New Zealand Bone and Mineral Society; September 2005; Perth, Australia.
14. Da Silva A.M.H., Alves J.M., Da Silva O.L., Da Silva Junior N.F. Two and Three-Dimensional Morphometric Analysis of Trabecular Bone Using X-Ray Microtomography (pCT). Revista Brasileira de Engenharia Biomedica, 2014, vol. 30, no. 2, pp. 93-101. DOI: 10.1590/rbeb.2014.011
15. Connor D.M., Hallen H.D., Lalush D.S., Sumner D.R., Zhong Z. Comparison of Diffraction-Enhanced Computed Tomography and Monochromatic Synchrotron Radiation Computed Tomography of Human Trabecular Bone. Physics in Medicine and Biology, 2009, vol. 54, no. 20, pp. 6123-6133. DOI: 10.1088/0031-9155/54/20/006
16. Gefen A. Finite Element Modeling of the Microarchitecture of Cancellous Bone: Techniques and Applications. Biomechanical Systems Technology, World Scientific Publishing, 2009, vol. 3, pp. 73-112. DOI: 10.1142/9789812771384_0003
17. Krestinina L.Y., Davis F.G., Schonfeld S, Preston D.L., Degteva M., Epifanova S., Akleyev A.V. Leukaemia Incidence in the Techa River Cohort: 1953-2007. British Journal of Cancer, 2013, vol. 109, pp. 2886-2893. DOI: 10.1038/bjc.2013.614
18. Degteva M.O., Tolstykh E.I., Vorobiova M.I., Shagina N.B., Shishkina E.A., Bougrov N.G. Techa River Dosimetry System: Present and Future. Voprosy radiacionnoy bezopasnosti, 2006, no. 1, pp. 81-95. (in Russian)
19. Volchkova A.Y., Chuvakova D.A., Shishkina E.A. Calculations of Tooth Enamel Doses from Internal Exposure Based on a Set of Voxel Phantoms by Example of the 1st Low Incisor. Radiat Safety Problems (Mayak Production Association Scientific Journal), 2009, no. 56, pp. 66-75. (in Russian)
20. Knuth D.E. The Art of Computer Programming, vol. II, Seminumerical Algorithms. N.Y., Addison-Wesley Longman, 1998.
Received March 12, 2018
УДК 004.942+001.851.573 Б01: 10.14529/ттр 180204
ПАРАМЕТРИЧЕСКАЯ СТОХАСТИЧЕСКАЯ МОДЕЛЬ
ГЕОМЕТРИИ КОСТИ
В.И. Заляпин1, Ю. С. Тимофеев2, Е. А. Шишкина2
1 Южно-Уральский государственный университет, г. Челябинск,
Российская Федерация
2
Российская Федерация
Целью настоящего исследования является разработка алгоритма параметрического моделирования костей с учетом особенностей их микроархитектуры, который позволяет генерировать фантомы гематопоэтических сегментов кости на основе литературных данных о микро- и макро- размерах сегментов. Мы предлагаем подход,
позволяющий моделировать кости путем их разбиения на небольшие сегменты, описываемые геометрическими фигурами простой формы, заполненные стохастически генерируемой стержнеподобной моделью трабекулярной структуры, с подходящим разрешением воксела. Предлагаемый метод позволяет избежать недостатков непараметрического моделирования, основанного на индивидуальной компьютерной томографии. Параметрический подход позволяет, также, моделировать индивидуальную изменчивость костно-специфических размеров. Модель, представленная в этой статье, может быть использована для описания геометрии гемопоэтических сайтов, которые, в свою очередь, служат основой для расчета доз облучения гемопоэтических клеток костного мозга от инкорпорированных бета-излучателей.
Ключевые слова: микро- и макро- структура трабекулярной кости; стохастическое моделирование; вокселизация.
Литература
1. Сох, L.G.E. Analysis of Bone Architecture Sensitivity for Changes in Mechanical Loading, Cellular Activity, Mechanotransduction, and Tissue Properties / L.G.E. Cox, B. van Rietbergen, C.C. van Donkelaar, K. Ito // Biomechanics and Modeling in Mechanobiology. -2011. - V. 10, № 2. - P. 701-712.
2. Vaughan, T.J. Multiscale Modeling of Trabecular Bone Marrow: Understanding the Micromechanical Environment of Mesenchymal Stem Cells During Osteoporosis / T.J. Vaughan, M. Voisin, G.L. Niebur, L.M. McNamara // Journal of Biomechanical Engineering. - 2015. - V. 137, № 1. - P. 011003.
3. McParland, B.J. Nuclear Medicine Radiation Dosimetry: Advanced Theoretical Principles /
B.J. McParland. - London: Springer Science and Business Media, 2010.
4. Dant, J.T. Monte Carlo Simulation of Age-Dependent Radiation Dose from Alpha- and Beta-Emitting Radionuclides to Critical Trabecular Bone and Bone Marrow Targets / J.T. Dant, R.B. Richardson, L.H. Nie // Physics in Medicine and Biology. - 2013. - V. 58, № 10. -P. 3301-3319.
5. Bolch, W.E. Hybrid Computational Phantoms for Medical Dose Reconstruction / W.E. Bolch,
C. Lee, M. Wayson, P. Johnson // Radiation Environmental Biophysics. - 2010. - V. 49, № 2. - P. 155-168.
6. Pafundi, D. An Image-Based Skeletal Dosimetry Model for the ICRP Reference Newborn -Internal Electron Sources / D. Pafundi, D. Rajon, D. Jokisch, C. Lee, W. Bolch // Physics in Medicine and Biology. - 2010. - V. 55, № 7. - P. 1785-1814.
7. Hough, M. An Image-Based Skeletal Dosimetry Model for the ICRP Reference Adult Male -Internal Electron Sources / M. Hough, P. Johnson, D. Rajon, D. Jokisch, C. Lee, W. Bolch // Physics in Medicine and Biology. - 2010. - V. 56, № 8. - P. 2309-2346.
8. O'Reilly, Sh.E. An Image-Based Skeletal Dosimetry Model for the ICRP Reference Adult Female - Internal Electron Sources / Sh.E. O'Reilly, L.S. DeWeese, M.R. Maynard, D. Rajon, M.B. Wayson, E.L. Marshall, WT. Bolch // Physics in Medicine and Biology. - 2016. - V. 61, № 24. - P. 8794-8824.
9. Fazzalari, N. Antero-Postero Differences in Cortical Thickness and Cortical Porosity of T12 to L5 Vertebral Bodies / N.L. Fazzalari, I.H. Parkinson, Q.A. Fogg, P. Sutton-Smith // Joint Bone Spine. - 2006. - V. 73, № 3. - P. 293-297.
10. Волчкова, А.Ю. Влияние формы кости на мощность дозы в костном мозге человека / А.Ю. Волчкова, Е.А. Шишкина, Е.И. Толстых// Адаптация биологических систем к природным и экстремальным факторам окружающей среды: Материалы VI Международной научно-практической конференции, Челябинск, ноябрь 8-9, 2016. - Челябинск, 2016. - С. 36-41._
Kg Bulletin of the South Ural State University. Ser. Mathematical Modelling, Programming
& Computer Software (Bulletin SUSU MMCS), 2018, vol. 11, no. 2, pp. 44-57
11. Roberts, S. On Some Geometric Properties of Human Ribs / S. Roberts, P. Chen // Journal of Biomechanics. - 1970. - V. 3. - P. 527-545.
12. Li, Z. Influence of Mesh Density, Cortical Thickness and Material Properties on Human Rib Fracture Prediction / Z. Li, M.W. Kindig, D. Subit, R.W. Kent // Medical Engineering and Physics. - 2010. - V. 32, № 9. - P. 998-1008.
13. Sutton-Smith, P. Quantitative Backscattered Electron Microscopy of Normal and Osteoporotic Bone / P. Sutton-Smith, H. Beard, N.L. Fazzalari // Proceedings of the Annual Scientific Meeting of the Australian and New Zealand Bone and Mineral Society. - Perth, 2005.
14. Da Silva, A.M.H. Two and Three-Dimensional Morphometric Analysis of Trabecular Bone Using X-ray Microtomography OCT) / A.M.H. Da Silva, J.M. Alves, O.L. Da Silva, N.F. Da Silva Junior // Revista Brasileira de Engenharia Biomedica. - 2014. - V. 30, № 2. - P. 93-101.
15. Connor, D.M. Comparison of Diffraction-Enhanced Computed Tomography and Monochromatic Synchrotron Radiation Computed Tomography of Human Trabecular Bone / D.M. Connor, H.D. Hallen, D.S. Lalush, D.R. Sumner, Z. Zhong // Physics in Medicine and Biology. - 2009. - V. 54, № 20. - P. 6123-6133.
16. Gefen, A. Finite Element Modeling of the Microarchitecture of Cancellous Bone: Techniques And Applications / A. Gefen // Biomechanical Systems Technology. - World Scientific Publishing, 2009. - V. 3. - P. 73-112.
17. Krestinina, L.Y. Leukaemia Incidence in the Techa River Cohort: 1953-2007 / L.Y. Krestinina, F.G. Davis, S. Schonfeld, D.L. Preston, M. Degteva, S. Epifanova, A.V. Akleyev // British Journal of Cancer. - 2013. - V. 109. - P. 2886-2893.
18. Дегтева, M.O. Дозиметрическая система бассейна реки Теча: настоящее и будущее / М.О. Дегтева, Е.И. Толстых, М.И. Воробьева, Н.Б. Шагина, Е.А. Шишкина, ПЛ . Бугров // Вопросы радиационной безопасности. - 2006. - № 1. - С. 81-95.
19. Волчкова, А.Ю. Расчет дозы внутреннего облучения зубной эмали на базе набора вексельных фантомов на примере 1-го нижнего резца / А.Ю. Волчкова, Д.А. Чувакова, Е.А. Шишкина // Проблемы радиационной безопасности (Научный журнал производственного объединения «Маяк>). - 2009. - № 56. - С. 66-75.
20. Knuth, D.E. The Art of Computer Programming. Vol. II, Seminumerical Algorithms / D.E. Knuth. - N.Y.: Addison-Wesley Longman, 1998.
Владимир Ильич Заляпин, К&НДИДсХТ физико-математических наук, профессор, кафедра «Математического анализа и методики преподавания математики>, ЮжноУральский государственный университет (г. Челябинск, Российская Федерация), zaliapinvi@susu.ru.
Юрий Сергеевич Тимофеев, кандидат технических наук, научный сотрудник, лаборатория биофизики, Уральский научно-практический центр радиационной медицины (г. Челябинск, Российская Федерация), swimmer_86@mail.ru.
Анатольевна Шишкина, К&НДИДсХТ биологических наук, старший научный сотрудник, лаборатория биофизики, Уральский научно-практический центр радиационной медицины (г. Челябинск, Российская Федерация), lena@urcrm.ru.
Поступила в редакцию 12 марта 2018 г.
