A graph model to describe the network connectivity of trabecular plates and rods

Introduction: The trabecular network is perceived as a collection of interconnected plate- (P) and rod-like (R) elements. Previous research has highlighted how these elements and their connectivity influence the mechanical properties of bone, yet further work is required to elucidate better the deeply interconnected nature of the trabecular network with distinct element formations conducting forces per their mechanical boundary conditions. Within this network, forces act through elements: a rod or plate with force applied to one end will transmit this force to a component connected to the other end, defining the boundary conditions for the loading of each element. To that end, this study has two aims: First, to investigate the connectivity of individually segmented elements of trabecular bone with respect to their local boundary conditions as defined by the surrounding trabecular network and linking them directly to the bone’s overall mechanical response during loading using a mathematical graph model of the plate and rod (PR) Network. Second, we use this model to quantify side artifacts, a known artifact when testing an excised specimen of trabecular bone, where vertical trabeculae lose their load-bearing capacity due to a loss of connectivity, ultimately resulting in a change of the trabecular network topology. Resuts: Connected elements derived from our model predicted apparent elastic modulus by fitting a linear regression (R 2 = 0.81). In comparison, prediction using conventional bone volume fraction results in a lower accuracy (R 2 = 0.72), demonstrating the ability of the PR Network to estimate compressive elastic modulus independent of specimen size or loading boundary condition. Discussion: PR Network models are a novel approach to describing connectivity within the trabecular network and incorporating mechanical boundary conditions within the morphological analysis, thus enabling the study of intrinsic material properties of trabecular bone. Ultimately, PR Network models may be an early predictor or provide further insights into osteo-degenerative diseases.

We will assume 26-connectivity between objects considered bone, called black points, and 6-connectivity between objects considered soft-tissue objects, called white points.
The following functions define distinct e-points and v-points based on their location which may be described uniquely by their three non-opposite neighboring s-points a, b, c of N p .
• e(a, b, p) = q|q ∈ N * p and 6-adjacent to both s-points a, b • v(a, b, c, p) = q|q ∈ N * p and 6-adjacent to e(a, b, p), e(b, c, p), and e(c, a, p) The following functions are used to expand the 3 * 3 * 3 neighborhood and can be defined in a 5 * 5 * 5 neighborhood of p.They are mainly used to preserve sharpness of corners during identification of surface voxels (Figure A.2).
• f 3 (a, b, p) = q|q / ∈ N p and 6-adjacent to f 2 (a, b, p) and f 2 (b, a, p).Above definitions for distinct points will be used frequently by the algorithm and do not change within N p and its expansion defined by f 1 , f 2 , f 3 .It is most efficient to compute them once and explicitly define them as a function of each distinct s-point combination.Once computed values can be saved in look-up tables and made easily accessible.
Thick-points: thick-points are part of a 2-voxel thick surface which can be identified using the definitions f 1 , f 2 and f 3 .Thick-points cannot be eroded during primary thinning, as peeling of a 2-voxel thick surface would result in a hole.

Surface points
The later described algorithm defines a peeling process where the surface of the 3D-Image is removed layer by layer.The surface of the image can be defined as the sum of the following s-open, e-open and v-open points: • A bone voxel is considered an s-open point if at least one s-point of N p is white before the iteration.

Shape preserving points
During primary thinning layers are removed until only shape preserving points remain and the 3D-image is reduced to an one-or two-voxel thick surface.In order d to identify whether a point is a shape preserving point in its neighborhood the following definitions may be used.
Within the 3 * 3 * 3 neighborhood a middle plane between two opposite s-points a and d can be defined in such way that it does not include a nor d.We will define the following conditions for all three pairs of opposite s-points (a, d), (b, e) and (c, f ).Let condition.
• C 1 be x ∈ {b, e, c, f }, • C 2 be x ∈ {b, e, c, f } and x, y are non-opposite, • C 4 be x ∈ {b, e, c, f }, x, y are non-opposite and x ∈ {p B , p S , p W }, • C 5 be x ∈ {b, e, c, f }, x, y are non-opposite and x, y ∈ {p B , p S , p W }.
A middle plane M(a, d, p) of N p in the 3 * 3 * 3 neighborhood and an extension of this plane EM(a, d, p) in the 5 * 5 * 5 neighborhood can be defined in the following way: We also define a surface of N p as follows: sur f ace(a, p) = {a} ∪ {e(a, x, p)|C 1 } ∪ {v(a, x, y, p)|C 2 }.
Shape preserving points fulfil the following conditions: • in N p two opposite s-points (a, d) exist such that EM(a, d, p) contains a 6-closed path of white points encircling p and each of sur f ace(a, p) and sur f ace(d, p) contains at least one black point before the iteration.
• in N p a pair of opposite s-points (a, d) exists such that d ∈ {p e , p s , p w }, a is white, e

Simple points
Erodible points are called simple points.The following definitions can be used first to identify simple points, second to classify the image later on.
Black objects ε(p), tunnels η(p) and cavities δ(p): We will use definitions of local topological parameters: the number of black objects ε(p), the number of tunnels η(p) and the number of cavities δ(p) within N p in order to classify distinct points of our 3D image (inner-surface points, outer-surface points, inner-arc points, arc-end points, surface-surface connection, surface-arc connection, arc-arc connection and isolated points).
• ε(p) is equal to the number of 26-connected black points within N * p • η(p) is one less than the number of 6-connected white components which intersect with an s-point of N p .
• δ(p) is always equal to 1 except when all s-points are black points.
The following four conditions were used to define a simple point: • K 2 be the number of tunnels in N * p is zeros (η(p)) • K 3 be N * p contains one 26-connected black component • K 4 be p has at least one black 26-neighbor

Classification of a 1-voxel thick skeleton
The following definitions will be used in order to classify the resulting 1-voxel thick skeleton and identify distinct geometric points.Commonly, the following types are distinguished: isolated points (I), inner-arc points (IA), edge-arc points (EA), arc-arc junctions (AA), inner-surface points (IS), edge-surface points (ES), surface-surface junctions (SS) and surface-arc junctions (SA).

Geometric classes:
As shown by Saha et al., these topological parameters (A.1.1)are independent of a surface within N p when the s-point within that surface is black.The remaining points will be called effective points.Each configuration of s-points results in a unique effective point configuration c e f f .The resulting configurations can be divided into 10 possible geometric classes depending on their s-point configuration.
Within each class, the geometric parameters are identical.These classes are mainly f used to help reducing computational cost.However, they do not directly respond to previously mentioned geometric points yet.
• Class 0: All s-points are bone.Number of effective voxels n e = 0.
• The number of possible effective point configurations is equal to 2 n e .The following parameters can be immediately identified for Class 0: ε

(p) η(p) δ(p), Class 1: ε(p) η(p) δ(p), Class 2: ε(p) η(p) δ(p).
For all other classes parameters can be computed once and saved within a lookup table for each class in order to speed up the code.

Subfields
The 3D-thinning algorithm describing a peeling process, where each iteration the outer layer is removed is inherently serial.However, it is possible to parallelize the point classification of erodible or shape preserving points, during each iteration.These points are evaluated based on their 26-neighborhood.Parallel removal of two g The abovementioned voxels are either eroded when they are considered simple points or preserved when they are considered a shape preserving point.At the beginning of the algorithm the value "0" is assigned to bone voxels and the largest negative 32-bit integer (−intmax) is assigned to marrow voxels.Zero-value voxels denote unmarked bone voxels, which are considered for erosion during each scan of the image.If an unmarked bone voxels is considered shape preserving it is marked, assigned the current iteration number i and cannot be deleted during primary thinning.Otherwise, if it is considered a simple point it is deleted by assigning the value thr = −maxint + i.As a result, one single physical image can be interpreted at different stages of the erosion process.Where simple points are always derived from the current image version I ≥ thr, shape preserving points are derived from the image at the beginning of the current iteration I > thr.The erosion process continues until two consecutive iterations contain the same voxel configuration, meaning no more voxel can be eroded or preserved.Primary thinning may results in a 2-voxel thick skeleton, which can be properly reduced using the following final thinning procedure.
Final thinning is a single iteration procedure reducing all thick points of a two voxel thick surface.During this final scan, all marked and unmarked bone voxel are considered for erosion.Thick points cannot be removed during primary thinning, as peeling of a two-voxel thick surface would remove the entire surface and may result in a false representation of the skeleton.

A.1.3. Classification
The resulting one voxel thick skeletons were classified using a three-dimensional digital topological characterization method [4].For each voxel, it is determined whether it is a surface voxel, a surface end voxel, an arc voxel, an arc-end voxel, an arc-arc intersection voxel, an arc-surface intersection voxel, a surface-surface intersection voxel, or an isolated voxel.A description of the basic principles of this algorithm can be found in A.1.1.
In accordance with Stauber et al. [7], slender planes were reduced to rods within two iterations and all rods shorter than four voxels were removed.Within this new optimized image individual rods and plates were separated by removing all points within a radius (r=2) of a junction point enabling individual labeling of rods with odd and plates with even indices.Using a 3D-bitmapping image dilation algorithm described by van den Boomgard et al. individually labeled images were remapped to their original geometry.
Within the original geometry, a point was identified as a junction point when at least one 26-adjacent voxel of one element was contained within a neighboring platej or rod-element resulting in a two-voxel thick junction surface between neighboring elements.Junctions were stored in a sparse, numeric adjacency matrix, where the location of each nonzero entry specifies a junction (edges) between two elements (nodes) building the basis of an unweighted network graph.

A.2. Spatial Network
A graph is a collection of nodes and edges with the following relationships: Nodes are vertices that correspond to objects and edges are connections between objects.Graph edges may have weights and may indicate the potential or strength of each connection.As these models are very general, the specific meaning of nodes, edges and weights highly depends on their application.Further, two types of graphs can be distinguished: undirected and directed graphs.In a directed graph, edges can only be traversed in one defined direction.Undirected graphs have edges that do not have a direction.The edges indicate a two-way relationship, in that each edge can be traversed in both directions.
For the spatial plate and rod network, nodes were defined as junction points between individual trabeculae and edges as trabeculae themselves.Therefore, weights can be chosen according to individual distinct properties (such as size, width, length, density) of each trabeculae.By applying common graph analysis algorithms such as the Boykov-Kolmogorov algorithm, which computes the maximum flow by constructing two search trees associated with nodes s and t, geometric parameters can be evaluated in regards to a specific loading case.In order to evaluate the flow through the whole specimen nodes s and t are artificial super-nodes connected to all nodes on the proximal and distal specimen surface respectively.
Figure A.1.: Illustration indicating the 26-neighbourhood of an exemplary point p.Each individual point is labeled according to their relative location top using cardinal points (N: North, E: East, S: South, W: West, T: Top, B: Bottom).Gray-scale code indicating 6 s-points, 12 e-points, and 8 vpoints making up the 26-neighborhood.

Class 1 :• Class 6 :• Class 7 :
Five s-points are bone.Number of effective voxels n e = 0. • Class 2: Two pairs of opposite s-points are bone.Number of effective voxels n e = 0. • Class 3: One pair of opposite and two non-opposite s-points are bone.Number of effective voxels n e = 1.• Class 4: One pair of opposite and another s-point are bone.Number of effective voxels n e = 2. • Class 5: Three non-opposite s-points are bone.Number of effective voxels n e = 4.One pair of opposite s-points are bone.Number of effective voxels n e = 4. Two non-opposite s-points are bone.Number of effective voxels n e = 7. • Class 8: Only one s-point is bone.Number of effective voxels n e = 12.• Class 9: No s-point is bone.Number of effective voxels n e = 20.

Table A
processing the entire traversal of the image.During each scan, the image border is identified as the aggregate of all s-open, e-open and v-open points.