Approximation of the Role of Mineralized Collagen Fibril Orientation in the Mechanical Properties of Bone: A Computational Study on Dehydrated Osteonal Lamellar Bone

Abstract

Bone is a natural composite of the hierarchical arrangement of mineralized collagen fibrils in various orientations. This study aims to understand how the orientation of the bone mineral, guiding the removal of water contained in the humidity-responsive layers during dehydration, affects its mechanical properties. A sublamellar pattern with mineralized collagen fibrils oriented between 0° to 150° at 5° angles was the model studied. Using basic transformational computational methods, dimensional change was calculated in the transverse and oblique planes of osteonal lamellar bone while considering bone components sensitive to dehydration in radial, tangential, and axial orientations. The anisotropy ratios of the change in the dimension of the variable mineralized collagen fibril orientations calculated using the computed model displayed values ranging between 0.847 to 2.092 for the transverse plane and 0.9856 to 1.0207 for the oblique plane. A comparison of the anisotropy results of the suggested model indicated that they approach the experimental results of both transversely and obliquely cut samples. As collagen fibril and mineral orientation take place both temporally and spatially in relationship with the static and dynamic loads placed on the different volumes of bone, the results may imply that the mechanical demands involved in bone resorption and deposition contribute to the formation of this multi-faceted and hierarchically structured natural composite.

1. Introduction

Bone is a natural composite material consisting of hierarchically structured mineralized collagen fiber. The organic and inorganic components of bone are contained within a hydrating medium and optimized to provide stiffness and toughness as a lightweight material [1,2,3,4,5,6,7,8]. Nature has organized these components hierarchically at different levels and length scales [9,10]. The organic component of bone, starting at the level of amino acids at the nanoscale, form the collagen fibrils [11,12], which are hierarchically organized further up as they associate with the inorganic component forming mineralized collagen bundles [9,13]. The micron-sized tessellated ellipsoidal mineral aggregates (made of about 10–100 mineral platelets of approximately 10 × 25 × 50 nm) [14,15] consist of laterally merging acicular crystals that are made of unit cells that fill in the gaps between the collagen fibrils [16,17]. These mineral crystals appear to align with the coiling fibril axis at about a 5° angle [18,19,20].

The anatomy of lamellar bone has been a long disputed topic with differing views on mineralized collagen fibril (MCF) organization as (a) uniform density rotated plywood-like MCF layers [21,22], (b) collagen-rich MCFs alternating with the collagen-poor MCF layers [23], (c) uniform density fibrils with apparently more or less concentrated mineralization [8,9], (d) collagen fibrils organized as both ordered unidirectional and fanning MCFs set at angles around the osteon, and less mineralized, disordered layers of MCFs that surround the irregular anatomical structures [9,13,24], and (e) as helical arrays of collagen fibrils appearing as twisted, curved, sinusoidal, coiled, oscillating and spiraling structures [25,26] to the extent that an osteon may be considered as concentric coils [9,13,24,27].

It is generally accepted that the collagen fibrils of the osteonal sublamellae rotate relative to the lamellar front, with the mineral planes [9,21] rotating relative to both the collagen fibrils as well as the lamellar front (Figure 1) [7,15,28]. Based on SEM images of mineralized bone samples, it is plausible to merge these approaches into an osteonal lamellar bone model that may explain the discrepancies between the various models. The model proposed here [29] consists of uniform density sublamellae of mineralized rotated collagen fibrils with their mineral planes oriented respectively at varying angles. The rotation of the MCFs displays a thin or thick appearance of the sublamellae as they are viewed in-plane or at an angle [3,29]. Therefore, such a lamellar pattern is visible under dry conditions, yet almost disappears under wet conditions [8].

The difference between the micrographic images of wet and dry bone have been investigated in terms of differences in dimension as well as mechanical properties in axial and radial directions [8,30,31,32,33]. Studies conducted on wet and dry MCFs display lateral contraction in the transverse plane with relatively less axial change in dimension [8,34,35,36]. This implies that the humidity-responsive layers, which are arranged parallel to the lamellar front, affect contraction [1,8,34]. In addition, the extent of mineralization and the orientation of mineral plates affect changes in dimension [8,29,37]. The organization of the humidity-responsive layers in bone contributes to bone anisotropy, indicating that the mechanical effect of water may not be isotropic [7,32,38,39] and that the orientation of the mineral plates within the hydrated collagen fibril may largely affect its mechanical properties [29,40,41].

This study examines the role of mineralized collagen fibril (MCF) orientation in influencing the mechanical properties of bone, particularly in dehydrated osteonal lamellar bone. The hypothesis is that the orientation of mineral crystals within MCFs guides the direction of removal of water from humidity-responsive layers, which significantly impacts bone’s mechanical behavior. A theoretical analysis of MCF orientations in transverse and oblique planes at the lamellar scale was conducted while considering the contribution of bone fluid to bone’s anisotropy and mechanical properties. Basic computational methods were used to approximate the dimensional changes of MCFs at 5° angles between 0° and 150°. These results were then compared with experimental data on changes in anisotropy between wet and dry osteonal lamellar bone.

2. Materials and Methods

2.1. Materials

A single lamella composed of sublamellar layers of cylindrical MCFs was used as the model. The MCFs were oriented at 5° angles between 0° and 150° (Figure 2). The (S-0) MCF was oriented along the z-axis [001] at 0° while the (S-90) MCF was oriented along the x-axis [100] at 90°. The (S-0) mineral plane (010) intercepted the y-axis, while the S-90 mineral plane (001) intercepted the z-axis (Figure 2). Each sublamella of a single, oriented MCF (numbered S-0, S-5, S-10, … through S-150) was layered along the y-axis on the x–z plane Figure 3a. The mineral plate was rotated at 5° intervals within the fibril of each sublamella along the x-, y-, and z- axes. This rotation pattern for the mineral ranged from 0° to 90° (S-0 to S-90) and was repeated at 5° intervals from 5° to 60° for the S-90 to S-150 collagen fibrils (Figure 3b) based on SEM images [7].

2.2. Methods

2.2.1. Transformation Matrix

Generalized average orientations and angular rotations for the MCFs were determined by transforming the fibril of the 0° sublamella (S-0) using transformation matrices and associated equations (Equations (1)–(8)) [40]. In the first stage of this study, the transformation process was conducted within the 0° to 150° range at 5° increments for $\theta$ to identify primary orientations and 3D transformations (expanded from the 120° range previously used by Utku (2021) [40]). Using this method, standard transformation matrix operations were applied to manipulate vectorial quantities. These operations are fundamental for translating, rotating, and scaling vectors in space and are widely used across various disciplines, such as mechanics, computer graphics, and robotics. The transformation matrices utilized here adhere to the principles of linear algebra, where matrices act as linear operators on vectors to transform them into new coordinate systems or configurations.

$$T_x\left(\theta \right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\theta & -sin\theta \\ 0& sin\theta & cos\theta \end{array}\right],$$

(1)

$$T_y\left(\theta \right)=\left[\begin{array}{ccc}cos\theta & 0& sin\theta \\ 0& 1& 0\\ -sin\theta & 0& cos\theta \end{array}\right],$$

(2)

$$T_z\left(\theta \right)=\left[\begin{array}{ccc}cos\theta & -sin\theta & 0\\ sin\theta & cos\theta & 0\\ 0& 0& 1\end{array}\right],$$

(3)

$$n_x=\left[100\right],$$

(4)

$$n_y=\left[010\right],$$

(5)

$$n_z=\left[001\right],$$

(6)

The orientation of MCFs was calculated as follows:

$$\mathrm{Collagen}\mathrm{orientation}=T_{\mathrm{y}}\ast n_{\mathrm{z}},$$

(7)

where T_x, T_y, and T_z denote the transformation matrices with respect to the x-, y-, and z-axes, while n_x, n_y, and n_z denote the vectors oriented along the x-, y-, and z-axes.

In the second stage of this study, the transformation process was repeated by rotating the initial angle of MCFs about the y-axis in 1° increments (from 0° to 10°, as S-0 + 1° to S-0 + 10°) through 150° (i.e., 150° + 10°) (Figure 2). This rotation (by ‘d’ degrees) was conducted based on findings indicating that the MCFs of the sublamellae closer to the Haversian Canal were not perfectly aligned vertically with the osteonal axis [9,13,24].

In the third stage of this study, in order to study the anisotropy of oblique planes, the orientations of the MCFs obtained from the second stage were rotated about the x-axis (by ‘f’ degrees) as depicted in Figure 2 and Figure 4a,b using Equations (9)–(11). The resultant oblique surface depicting the MCF orientations has been given in Figure 5a–c.

2.2.2. Calculation of Dimensional Change in Transverse and Oblique Planes

A transformation process was conducted to determine each mineral plate’s surface normal because it was assumed that as the organic component was dehydrated, collagen contracted towards the mineral plate. This study aimed to correlate the orientation of bone mineral (which guides the removal of water contained in the humidity-responsive layers during dehydration) with the anisotropy of its mechanical properties. Thus, the sublamellar pattern of mineralized collagen fibrils given in Section 2.2.1 was used to calculate the three-dimensional change in bone components sensitive to dehydration. Here, with the removal of water mainly from the organic collagen matrix, a dehydration-led contraction in collagen diameter was observed towards the mineral plate, i.e., in the direction normal to the mineral plate’s surface. Thus, the dehydration-led reduction in dimension was vectorial and named as the contraction vector [3,8,29]. The contraction direction was calculated according to Equation (8) and displayed in Figure 6.

$$\mathrm{Contraction}\mathrm{in}\mathrm{mineralized}\mathrm{collagen}\mathrm{fibril}=T_{\mathrm{y}}\ast T_{\mathrm{x}}\ast T_{\mathrm{z}}\ast n_{\mathrm{z}},$$

(8)

The axial projections of the contraction vector of the dehydrated fibril were calculated as u_t, v_t, and w_t for the x-, y-, and z-axes using the transformation matrices given in (Equations (1)–(8)). Since contraction is manifested as displacement, the absolute values of displacements at 5°, 10°, 15°, 20°, 25°, 30° or 40° angles (‘a’ degree) were summed up to 120°, 125°, 130°, 135°, 140° or 150° (‘b’ degrees). The contraction sums (U_at, V_at, and W_at) denote the axial sums of the contraction vector in the respective orthogonal axes.

To study the possible contraction in an oblique plane, the model was assumed to have rotated about the x-axis by ‘f’ degrees, ∼40° (−40° ± 5°) the angle described by Faingold et al. [41] (Figure 4b). The total contraction was calculated using Equations (9)–(11), and denoted as U_ao, V_ao, and W_ao, the total magnitude of the contraction vector in the respective global axes.

$${\mathrm{U}}_{\mathrm{ao}}=T_{\mathrm{x}}\ast n_{\mathrm{x}},$$

(9)

$${\mathrm{V}}_{\mathrm{ao}}=T_{\mathrm{x}}\ast n_{\mathrm{y}},$$

(10)

$${\mathrm{W}}_{\mathrm{ao}}=T_{\mathrm{x}}\ast n_{\mathrm{z}},$$

(11)

A summary of the methods and variables used in this study is listed in

Table 1. A summary of the methods and variables used in this study.

Stage of Study	Angle Title	Axial Rotation About	Angle Range	Step Size	Angles
1st	‘a’	y-axis	0° to 150°	5°	5°, 10°, 15°, 20°, 25°, 30°, 40°
	‘b’	y-axis	120° to 150°	5°	120°, 125°, 130°, 135°, 140°, 150°
2nd	‘d’	y-axis	1° to 10°	1°	1°, 2°, 3°, …, 10°
3rd	‘f’	x-axis	35°− 45°	1°	−35°, −36°,…, −45°

.

2.2.3. The Anisotropy Ratio

The anisotropy ratios (AnR) were calculated as the ratio of contraction (experimentally as microhardness or stiffness) of the transverse (cross-sectional) plane (obtained parallel to the osteonal axis (W)) to those obtained perpendicular to the osteonal axis (V) (Figure 4a).

The transverse AnR was calculated as follows:

$$An{R}_{\mathrm{ct}}=W_{\mathrm{at}}/V_{\mathrm{at}}=1/An{R}_{\mathrm{rt}}=1/(V_{\mathrm{at}}/W_{\mathrm{at}}),$$

(12)

while the oblique AnR_co was calculated as follows:

$$An{R}_{\mathrm{co}}=W_{\mathrm{ao}}^{\prime }/V_{\mathrm{ao}}^{\prime }=1/An{R}_{\mathrm{ro}}=1/(V_{\mathrm{ao}}^{\prime }/W_{\mathrm{ao}}^{\prime }),$$

(13)

where a, t, o, c, and r denote the angle, transverse plane, oblique plane, compliance, and resistance, respectively. AnR_ct and AnR_co denote anisotropy due to the compliant element (c) in transverse and oblique planes, while the AnR_rt and AnR_ro denote anisotropy due to the contraction resistive (r) element in transverse and oblique planes. The AnR_ct was considered as the change in dimension of the organic component of MCF that is associated with water. Thus, contraction due to dehydration involved the compliant element (c) as the primary property of wet bone, i.e., as the capacity of the fibrils to contract or deform. On the contrary, the AnR obtained with dehydration was described as the AnR of contraction-resistive (r) element (AnR_rt), being the inverse of the compliant element (AnR_ct) [29].

2.3. Statistical Analysis

The transverse and oblique plane results were analyzed using an independent t-test for two data sets with equal variance (homoscedastic) for statistical variability at a significance level of 0.05. For the transversely cut bone, the AnR_ct results of each ‘a’ angle up to ‘b’ degrees were compared through ‘d’ degrees. As an example, 5° angles up to 120° degrees were compared with 5° angles up to 130° from 0° through 10°. For the obliquely cut bone, the AnR_co results of each ‘a’ angle up to ‘b’ degrees were compared through ‘f’ degrees. As an example, 5° angles up to 120° degrees were compared with 5° angles up to 130° degrees from 35° through 45°.

3. Results

The dimensional change due to the contraction of MCFs was studied at the lamellar length scale of an osteonal bone using a model of MCFs angulated at 5° angles from 0° up to 150°. The projections of the vectors in the x-, y-, and z-axes were evaluated at 5°, 10°, 15°, 20°, 25°, 30° or 40° intervals (‘a’ angle). The MCF orientations and the contraction vector projections with respect to global coordinates were calculated (

Table 2. Axial intercepts of MCF orientations and the contraction unit vectors per sublamellar layer.

	Collagen Fibril Orientation			Contraction Vector Projections
Sublamellae	X-Axis	Y-Axis	Z-Axis	u	v	w
S-0	0.0000	0	1.0000	0.0000	1.0000	0.0000
S-5	0.0872	0	0.9962	0.0793	0.9924	0.0941
S-10	0.1737	0	0.9848	0.1413	0.9699	0.1986
S-15	0.2588	0	0.9659	0.1853	0.9330	0.3085
S-20	0.3420	0	0.9397	0.2115	0.8830	0.4190
S-25	0.4226	0	0.9063	0.2212	0.8214	0.5257
S-30	0.5000	0	0.8660	0.2165	0.7500	0.6250
S-35	0.5736	0	0.8192	0.2004	0.6710	0.7138
S-40	0.6428	0	0.7660	0.1759	0.5868	0.7904
S-45	0.7071	0	0.7071	0.1465	0.5000	0.8536
S-50	0.7660	0	0.6428	0.1152	0.4132	0.9033
S-55	0.8192	0	0.5736	0.0850	0.3290	0.9405
S-60	0.8660	0	0.5000	0.0580	0.2500	0.9665
S-65	0.9063	0	0.4226	0.0359	0.1786	0.9833
S-70	0.9397	0	0.3420	0.0194	0.1169	0.9930
S-75	0.9659	0	0.2588	0.0085	0.0670	0.9977
S-80	0.9848	0	0.1737	0.0026	0.0302	0.9995
S-85	0.9962	0	0.0872	0.0003	0.0076	0.9999
S-90	1.000	0	0.0000	0.0000	0.0000	1.0000
S-95	0.9962	0	0.0872	0.0793	0.9924	0.0941
S-100	0.9848	0	0.1737	0.1413	0.9699	0.1986
S-105	0.9659	0	0.2588	0.1853	0.9330	0.3085
S-110	0.9397	0	0.3420	0.2115	0.8830	0.4190
S-115	0.9063	0	0.4226	0.2212	0.8214	0.5257
S-120	0.8660	0	0.5000	0.2165	0.7500	0.6250
S-125	0.8192	0	0.5736	0.2004	0.6710	0.7139
S-130	0.7660	0	0.6428	0.1759	0.5868	0.7904
S-135	0.7071	0	0.7071	0.1465	0.5000	0.8536
S-140	0.6428	0	0.7660	0.1152	0.4132	0.9033
S-145	0.5736	0	0.8192	0.0850	0.3290	0.9405
S-150	0.5000	0	0.8660	0.0580	0.2500	0.9665

). The axial projections of the contraction vectors (u_a, v_a, and w_a) were added to give the total axial contraction (U_at, V_at, W_at) up to ‘b’ angle for the transverse plane (

Table 3. Transverse plane MCF contractions summed up to 120°–150° with 0° to 10° rotation about the y-axis.

		120°	125°	130°	135°	140°	150°
ua		0.2165	0.2004	0.1759	0.1465	0.1152	0.0580
va		0.7500	0.6710	0.5868	0.5000	0.4132	0.2500
wa		0.6250	0.7139	0.7904	0.8536	0.9033	0.9665
∑5°	Ub	2.9576	3.1579	3.3338	3.4803	3.5955	3.7385
∑5°	Vb	14.8497	15.5207	16.1075	16.6075	17.0207	17.5997
∑5°	Wb	15.4833	16.1972	16.9876	17.8411	18.7445	20.6515
∑10°	Ub	1.5097	-	1.6856	-	1.8008	1.8588
∑10°	Vb	7.6029	-	8.1897	-	8.6029	8.8529
∑10°	Wb	8.1378	-	8.9282	-	9.8316	10.7981
∑15°	Ub	1.0166	-	-	1.1631	-	1.2211
∑15°	Vb	5.1830	-	-	5.6830	-	5.9330
∑15	Wb	5.6847	-	-	6.5383	-	7.5048
∑20°	Ub	0.8058	-	-	-	0.9210	-
∑20°	Vb	4.4698	-	-	-	4.8830	-
∑20°	Wb	3.9990	-	-	-	4.9023	-
∑25°	Ub	-	0.6865	-	-	-	0.7445
∑25°	Vb	-	3.9424	-	-	-	4.1924
∑25°	Wb	-	3.3392	-	-	-	4.3058
∑30°	Ub	0.4910	-	-	-	-	0.5490
∑30°	Vb	2.7500	-	-	-	-	2.0000
∑30°	Wb	3.2165	-	-	-	-	4.1830
∑40°	Ub	0.3950	-	-	-	-	-
∑40°	Vb	2.3670	-	-	-	-	-
∑40°	Wb	2.4149	-	-	-	-	-

). The AnR_ct for the transverse plane (

Table 4. The AnR_ct for the Transverse Plane Rotated by 0° to 10° (‘d’ angle) about the y-axis.

‘a’	‘b’	0°	1°	2°	3°	4°	5°	6°	7°	8°	9°	10°
∑5°	120°	1.043	1.084	1.127	1.171	1.216	1.264	1.312	1.363	1.415	1.470	1.526
∑5°	130°	1.055	1.097	1.141	1.186	1.233	1.282	1.333	1.385	1.439	1.495	1.554
∑5°	140°	1.101	1.146	1.192	1.239	1.289	1.340	1.393	1.449	1.506	1.565	1.627
∑5°	150°	1.173	1.220	1.269	1.320	1.373	1.427	1.484	1.543	1.604	1.668	1.734
∑10°	120°	1.070	1.112	1.155	1.200	1.246	1.293	1.342	1.393	1.446	1.500	1.557
∑10°	130°	1.090	1.133	1.178	1.224	1.272	1.321	1.372	1.425	1.480	1.537	1.596
∑10°	140°	1.143	1.188	1.235	1.284	1.334	1.386	1.440	1.496	1.555	1.615	1.678
∑10°	150°	1.220	1.268	1.318	1.369	1.423	1.479	1.537	1.597	1.659	1.724	1.791
∑15°	120°	1.097	1.139	1.182	1.227	1.273	1.320	1.370	1.421	1.473	1.528	1.584
∑15°	135°	1.150	1.195	1.241	1.289	1.339	1.390	1.443	1.498	1.555	1.614	1.675
∑15°	150°	1.265	1.314	1.365	1.417	1.472	1.528	1.587	1.648	1.711	1.776	1.844
∑20°	120°	0.895	0.931	0.969	1.009	1.049	1.091	1.135	1.180	1.227	1.275	1.325
∑20°	140°	1.004	1.045	1.087	1.130	1.176	1.223	1.271	1.322	1.374	1.429	1.485
∑25°	125°	0.847	0.883	0.919	0.957	0.997	1.038	1.080	1.124	1.169	1.216	1.265
∑25°	150°	1.027	1.068	1.110	1.155	1.200	1.248	1.297	1.348	1.401	1.456	1.514
∑30°	120°	1.170	1.211	1.254	1.298	1.344	1.391	1.439	1.489	1.540	1.592	1.647
∑30°	150°	2.092	2.180	2.272	2.368	2.469	2.575	2.685	2.800	2.920	3.045	3.176
∑40°	120°	1.020	1.059	1.098	1.140	1.182	1.226	1.271	1.318	1.367	1.417	1.468

) and AnR_co for the oblique plane were calculated (

Table 5. The AnR_co values (W/V) for the oblique plane for the y-axis rotation of 0° (x-axis rotation by −35° to −45° (‘f’ angle)).

‘a’	‘b’	−35°	−36°	−37°	−38°	−39°	−40°	−41°	−42°	−43°	−44°	−45°
∑5°	120°	1.0074	1.0066	1.0059	1.0051	1.0044	1.0037	1.0029	1.0022	1.0015	1.0007	1.0000
∑5°	130°	1.0094	1.0085	1.0075	1.0066	1.0056	1.0047	1.0037	1.0028	1.0019	1.0009	1.0000
∑5°	140°	1.0171	1.0154	1.0136	1.0119	1.0102	1.0085	1.0068	1.0051	1.0034	1.0017	1.0000
∑5°	150°	1.0285	1.0256	1.0227	1.0198	1.0169	1.0141	1.0112	1.0084	1.0056	1.0028	1.0000
∑10°	120°	1.0121	1.0108	1.0096	1.0084	1.0072	1.0060	1.0048	1.0036	1.0024	1.0012	1.0000
∑10°	130°	1.0153	1.0138	1.0122	1.0107	1.0091	1.0076	1.0061	1.0045	1.0030	1.0015	1.0000
∑10°	140°	1.0238	1.0213	1.0189	1.0165	1.0141	1.0117	1.0094	1.0070	1.0047	1.0023	1.0000
∑10°	150°	1.0355	1.0319	1.0282	1.0246	1.0210	1.0175	1.0139	1.0104	1.0069	1.0035	1.0000
∑15°	120°	1.0164	1.0147	1.0131	1.0114	1.0098	1.0081	1.0065	1.0049	1.0032	1.0016	1.0000
∑15°	135°	1.0250	1.0224	1.0199	1.0173	1.0148	1.0123	1.0098	1.0074	1.0049	1.0024	1.0000
∑15°	150°	1.0421	1.0378	1.0334	1.0291	1.0249	1.0207	1.0165	1.0123	1.0082	1.0041	1.0000
∑20°	120°	0.9806	0.9825	0.9845	0.9864	0.9884	0.9903	0.9923	0.9942	0.9961	0.9981	1.0000
∑20°	140°	1.0007	1.0006	1.0006	1.0005	1.0004	1.0003	1.0003	1.0002	1.0001	1.0001	1.0000
∑25°	125°	0.9712	0.9741	0.9770	0.9799	0.9827	0.9856	0.9885	0.9914	0.9942	0.9971	1.0000
∑25°	150°	1.0047	1.0042	1.0038	1.0033	1.0028	1.0023	1.0019	1.0014	1.0009	1.0005	1.0000
∑30°	120°	1.0280	1.0251	1.0222	1.0194	1.0166	1.0138	1.0110	1.0082	1.0055	1.0027	1.0000
∑30°	150°	1.1328	1.1185	1.1044	1.0906	1.0771	1.0637	1.0506	1.0377	1.0250	1.0124	1.0000
∑40°	120°	1.0035	1.0032	1.0028	1.0025	1.0021	1.0018	1.0014	1.0011	1.0007	1.0004	1.0000

Calculations have been provided only for the 0° rotation angle.

).

3.1. Transverse Plane

The axial intercepts of the collagen fibril orientations and MCF contraction vectors starting with S-0 at 0° angle, at 5° intervals up to 150° in the transverse plane are listed in Table 2. The axial intercepts of the contraction vectors, summed at ‘a’ angle intervals up to ‘b’ angles ranging between 120° and 150° are listed in Table 3. In the transverse plane, U_at displayed minimum contraction at all ‘a’ angles. At 5°, 10°, 15°, 30° and 40° angles W_at, and at 20° and 25° angles V_at displayed maximum contraction. Transverse plane results displayed an inverse relationship between the results of the projections in the y- and z-axes. The smallest contraction being U_at, the results of V_at and W_at demonstrated the behavior of a positive Poisson’s ratio material, where radial expansion is expected to constrict or minimize axial displacement (as was shown in previous studies [8,29]). The radial–tangential vs. axial, or the radial–axial vs. tangential correlation displayed here may be complemented by the strains [42] and inherent pre-stresses resulting from the coiling MCFs within the respective sublamella [36,43]. The AnR_ct (W/V) for the transverse plane was calculated (Table 4) and compared with the experimental dry and wet bone stiffness AnR results [41].

In the second stage of this study, the axial intercepts of MCF contractions starting with (S-0) set at 0° angle were rotated at a step size of 1° up to a maximum of 10° (by ‘d’ degrees). The ‘a’ angle intervals were summed from 120° (121°…130°) up to 160° (151° …160°). The results of this stage are lengthy and are available in the Supplementary Materials. The transverse plane AnR results for the 0° to 10° angle range are listed in Table 4. At the y-axis rotation angle of 0°, the lowest AnR_ct was 0.847 for the MCF contraction summed up to 125° at 25° angles, while the highest AnR_ct was 3.176 for the MCF contraction summed up to 150° at 30° angles. There was a gradual increase in AnR_ct (a) as the rotation angle (‘d’) was increased from 0° to 10° and also (b) as the angle of summation (angle ‘b’) was increased from 120°+ to 150°+. The 20° angles summed up to 120° versus 140° indicated that the AnR_ct increased from 1.325 to 1.485. Theoretically, the increase in AnR_ct was due to a higher number of axially oriented contractions accompanied by a reduced radial (or lateral) contraction. In other words, the ratio of MCFs oriented along the x-axis (with their contraction vectors oriented along the z-axis) to those oriented along the z-axis (with their contraction vectors oriented along the x-axis) was increased. Experimentally attained hardness AnR of 0.80 for wet bones [41] indicated that the wet bone stiffness was lower in the axial direction than the radial. In this study, the 20° and 25° MCF orientations (summed up to 120°) correlated with the experimental AnR < 1 for wet bone [41]. The experimental AnR < 1 was attained for the 20° MCF orientations (summed up to 120°) approximately between 0°–3° (‘d’) degrees, and for the 25° MCF orientations (summed up to 120°) between 0°–5° (‘d’) degrees.

3.2. Oblique Plane

Oblique plane MCF contractions were derived from the transverse plane results (given in Table 4) and listed in Table 5. The experimental values had been obtained at an oblique angle of ∼40° [41]; therefore, the transverse plane was rotated about the x-axis by −35° to −45° (‘f’) at a step size of −1°. Table 5 lists the oblique plane AnR_co results for only the 0° y-axis rotation (‘d’). The complete AnR_co data for the 0° to 10° (‘d’) is provided in the Supplementary Materials. For this angle range, the computed AnR_co varied between 1.0000 and 1.0421. For the 20° angle (summed up to 120°) and 25° angle (summed up to 125°), the AnR_co increased from −35° to −45°. However, for the other ‘a’ angles, the AnR_co decreased within that angle range.

The AnR_ct and AnR_co data were statistically evaluated using an independent t-test to determine the significance of the difference between the AnR results as given in Section 2.3 (

Table 6. T-test results (p_ct/p_co) for the transverse plane AnR_ct (0° through 10°) versus the oblique plane AnR_co (−35° through −45°) (p stands for probability).

Angle a	Angle b	120°	125°	130°	135°	140°
∑5°	130°	0.3935/0.2025
∑5°	140°	0.1373/0.0085		0.2132/0.0319
∑5°	150°	0.0183/0.0010		0.0323/0.0025		0.1328/0.0529
∑10°	130°	0.3442/0.2095
∑10°	140°	0.1283/0.0211		0.1902/0.0778
∑10°	150°	0.0104/0.0029		0.0246/0.0090		0.1224/0.0952
∑15°	135°	0.1683/0.0863
∑15°	150°	0.0585/0.0055			0.0443/0.0499
∑20°	140°	0.0265/0.00002
∑25°	150°		0.0018/0.00001
∑30°	150°	0.00005/0.0006

T-test results for the AnR_co have been provided only for the 0° (‘d’) rotation angle. Statistically insignificant values have been shown in bold.

). In general, the t-test results indicated that the data sets were not significantly different at smaller ‘a’ angles (namely, 5°, 10°, and 15°) when the ‘b’ summation angles were closer to each other (namely, 120° vs. 130° than 120° vs. 150°). The transverse plane AnR_ct results that compared ‘d’ rotation about the y-axis indicated that there was lower significant difference between the results of 120°, 130°, and 140° (‘b’) summation angles. The oblique plane AnR_co results that compared ‘f’ rotation about the x-axis indicated that the results were not significantly different between 120° and 130° (‘b’) summation angles. However, at higher ‘a’ angles (20°, 25°, 30°) the differences between 120° and 130°–150° angles were statistically significant.

4. Discussion

In this study, basic structural elements of bone have been modeled using specific boundary conditions in order to understand the possible organization of its components. To this aim, the mechanical changes that the osteonal bone undergoes through dehydration–rehydration were investigated using a 3D transformation code that determined anisotropy [29]. While a 5° angular interval was maintained, the angular range (‘b’) was expanded from 120° to 150° [8,29] and complemented with a study on oblique plane anisotropy. The results were evaluated in light of the experimental findings [41] on bone anisotropy at the fibrillar, lamellar, and osteonal length scales [7,32,36,39,44].

In a previous study, the computed AnR_ct (W/V) of 0.89 and AnR_rt (V/W) 1.12 were obtained using the axial projections of MCFs at 20° angle intervals [29]. These ratios appeared to correlate with the experimental findings at the lamellar length scale (stiffness AnR for wet and dry bone being 0.80 and 1.13, respectively) [41]. The experimental stiffness AnR < 1 for wet bone (axial/radial) meant that the axial stiffness was lower than the radial stiffness. As such, it indicated either (i) the presence of substantial axial (W) arrangement of contractible elements of MCF, making wet bone more compliant in that direction, or (ii) increased radial resistance exerted by both the hydraulic strength of bone fluid and mineral orientation [45].

The transverse plane AnR_ct results indicated that at all ‘a’ angles (5°, 10°, 15°, 20°, 25°, 30°, and 40°), summations (‘b’) up to 120° displayed ratios that were closer to the values stated in the literature. At higher summation (‘b’) angles (130°, 140°, and 150°) as well as higher y-axis rotation angles (‘d’), the AnR_ct were higher than the experimental results. The computed AnR_ct values as high as 3.176 (obtained at 30° summed up to 150°) demonstrated the extent of the effect of MCF orientation and organization on mechanical properties. Within the results listed in Table 4, only the AnR_ct for the 20° and 25° angles (summed up to 120° and 125° respectively) displayed values closer to that of wet osteonal lamellar bone (AnR = 0.80). The fact that these values were obtained at 0°–5° (‘d’) rotation appears to correlate with the experimental findings [9,13,24], which state that the MCF orientation may not be perfectly aligned with the osteonal axis at 0°. Experimental studies on the oblique-cut ($\theta$ ≈ 40°) osteonal bone samples indicated that stiffness AnR did not vary significantly between dry and wet samples. An AnR of 1 for these bones implied planar isotropy across the oblique plane. In this study, similar oblique plane results were obtained computationally by rotating the transverse plane about the x-axis (‘f’).

Statistical evaluations of the data using an independent t-test indicated that the data sets were not significantly different at smaller ‘a’ angles (namely, 5°, 10°, and 15°) when the ‘b’ summation angles were closer to each other (namely, 120° vs. 130° than 120° vs. 150°). The AnR_ct t-test results (that compared ‘d’ rotation data) indicated that the results were not significantly different between 120°, 130°, and 140° (‘b’) summation angles. The AnR_co t-test results (that compared ‘f’ rotation data) indicated that the results were not significantly different between 120° and 130° (‘b’) summation angles. However, for the 20°, 25°, and 30° angles, the differences between the 120° vs. 130°–150° summation angles were statistically significant. In an axially loaded bone, in order to ensure bone integrity, a specific AnR_ct range must be maintained despite the change in the orientation of MCFs. The AnR_co significance, observed in the oblique plane only between 120° and 130° summation angles, implies that the MCFs may have a narrow AnR_co range, thus making bone more “liable” to fractures in that orientation. The oblique plane, in fact, correlates with the surface exposed to shear stresses, which leads to the commonly observed shear fractures.

Experimentally and computationally, MCF orientation has been shown to be significant in the determination of the mechanical properties of bone. In previous studies, using computational methods and models with different MCF orientations, AnR_ct was shown to vary from 0.1527 to 6.5505, pointing to the significance of MCF orientation in determining the mechanical properties and bone anisotropy even in the very limited length scale considered here [8,10,29,40]. Using the change in mechanical properties of bone with dehydration, the stiffness AnR has been shown to vary at different length scales from fibril to lamellae to osteon, with higher stiffness in dry samples than in wet ones on all length scales [41]. The computed results indicate that angulations up to 120° give AnR_ct that correlate with the experimental AnR results. An axial to transverse strength ratio of 2.7:1 for tension and 1.3:1 for compression in mechanically tested lamellar bone samples [46] or the experimental AnR of 1.5, obtained on the length scale of the whole bone [47], have been approached computationally at higher summation (‘b’) and rotation angles (‘d’) (Table 3 and Table 4). The highest (and the most unrealistic) computational AnR_ct of 3.176, attained at 30° summed up to 150°, implied that MCF orientation could be an important indicator of the mechanical integrity of bone. Specifically, the results of the model, which have not been experimentally verified, pointed to the possible range of MCF orientations.

One of the possible complications affecting the significance of experimental results might be the difficulty of obtaining samples cut at exactly the same transverse or oblique angle during sample preparation. This is specifically important as the fibril orientation within the osteon varies from axially oriented (the lamellae closer to the Haversian Canal) to angulated (the lamellae closer to the cement line) [24]. Therefore, the mineral plates of the dehydrated collagen fibrils, which are layered at various angles ranging from 0° to 150° (through S-0 to S-150), may not only determine the contraction vector but also act as local constrictors within the tissue [45]. The sampling length [48], pre-stresses generated in bone [36] and tissue response due to the axial (W_at), tangential (U_at) and radial displacements (V_at) also emerge as critically important elements which determine the mechanical properties of this multi-faceted and hierarchically structured material. Using biomimicry, the design of bio-inspired materials with specific ‘a’ and ‘d’ angulation and ‘b’ summation of the MCFs may be fine-tuned to fulfill the necessary structure-location demands. Based on strain-dependent properties, bone structure may be statically and dynamically “programmed” to remain within the mechanical properties that ensure bone integrity. Bone responds to axial mechanical overloads using anisotropic deformation mechanisms, forming kinks and shear cracks at the interface of ordered and disordered regions, and to radial loads by forming cracks in order to reduce the local fracture resistance [46]. Within the concentric lamellar construction of the osteon, the structural design may be further complicated by both the osteon size [44] and the variation in MCF orientation across the osteon [24].

The organization of fibrils in bone takes place both temporally and spatially in response to the static and dynamic loads placed on the different sections of bone. Local mechanical strain patterns guide the deposition of collagen fibrils in various orientations that reinforce the tissue against local deflections [49,50]. The deposition and organization of the collagen bundles are statically determined as they are extruded from the cell. Non-uniform mineralization of the deposited collagen fibrils may also be observed [51]. Therefore, fibril orientation and the thickness of sublamaellae may be more or less constant through the lamellae or the osteonal volume due to the variations in mechanical strains observed at different voxels of the bone sample. This makes bone a strain/modulus location-dependent composite material with the advantage of structural modification based on local demands of anisotropy and transverse isotropy [52].

Another aspect that may be considered with this model is the location of the mineral: intra-fibrillar versus extra-fibrillar. According to the Hodge–Petruska model, considering that the 60% of each axial period is a gap region, the average volume fraction of intra-fibrillar mineral may be 0.43 with an upper possible limit of 0.56 in fully mineralized cortical bone [12]. Although the model takes dominantly the intra-fibrillar mineral, based on SEM images, we consider that the extra-fibrillar mineral tessellations may follow the dominant local direction of contraction of the MCFs [14]. Given that the individual slightly elliptical cross-sectioned collagen fibrils are bridged within a close-packed configuration, they have been thought to act as 3D volume fillers [14]. Therefore, the contribution of extra-fibrillar mineral platelets may be considered in two ways: (i) in terms of the collagen fibril orientation deposited according to the local strain distribution, (ii) the differences in the mechanical properties of the mineral platelets in the various azimuthal rotations [45,53].

In addition to the mechanically induced changes, chemically induced changes are important as they may affect bone integrity in diseased and aged bone [46]. Up to now, these computations have been used to study the mechanical consequences of chemically induced changes in terms of anisotropy in dry versus wet bone samples [41]. Using the mechanical consequences, both the orientation and the frequency of the MCFs have been approximated, implying the presence of a structure-location-modulus relationship within the bone. Thus, the stiffness (and hardness) AnR results of a sample may enable us to estimate its probable structural organization as well as its location.

Our model aimed to study a single lamella, which has been described as a rotated plywood structure (with uniform collagen density) defined by specific angles of tilt between lamellar sublayers, with a characteristic “back-flip” and different azimuthal rotations of the mineral platelets [26]. Therefore, imitating the back-flip sublamella, we summed up to at least 120° angle and further continued to sum up to 150°. Our model attempted to understand the mechanical properties of bone displayed in the literature [8,41]. In that study [8], the lamellae located in the central region of an osteon, away from both the cementing line and the Haversian canal, were sampled. Therefore, it would be less likely to observe oscillating plywood structures, which may be observed closer to the Haversian canal [26]. However, the model can be used to study both the oscillating plywood (MCFs oscillating within 30°) and twisted plywood (MCFs rotating within 180°) structures. The oscillating plywood structures would be expected to have a lower AnR_ct and thus a higher AnR_rt than the results obtained here. This would imply that the MCFs would be organized along the osteonal axis. On the other hand, twisted plywood structures would be expected to display higher AnR_ct and thus lower AnR_rt values.

In summary, at the lamellar length scale, the overall sublamellar patterning may be organized at ‘a’ angles up to ‘b’ summations and with ‘d’ rotations, which provide not only static but also hydraulic strength with its axially and radially arranged contracting elements in the sublamellae that contribute to the mechanical properties of osteonal lamellar bone. Our results appear to support the presence of MCF orientations at lower ’a’ angles (5° to 20°) summed up to lower ’b’ angles (120°–130°). This study is significant in the demonstration of the contribution of MCF orientation to bone mechanical properties.

5. Conclusions

The results of this study indicate the presence of tangentially, axially, and radially oriented bone components sensitive to hydration. Such sublamellar patterning provides both static and hydraulic strength to maintain bone’s structural and physical integrity, as it affects the anisotropy in contraction/expansion properties. A comparison of the anisotropy results of the suggested bone model approaches the experimental results both in transversely and obliquely cut samples. The results of this study demonstrate the structure–modulus–location relationship within bone. As has been established, mineralized collagen fibril (MCF) organization is a critical determinant of bone mechanics to the extent that maintaining anisotropy ratios that are within the experimentally determined limits may be a critical factor in bone integrity during bone resorption and deposition. MCF organization constructs bone’s hierarchical architecture while seeking adequate toughness and strength. This study demonstrates the importance of bone’s hierarchical microstructure and its effect on its mechanical properties as a load-bearing material (i.e., the structure-property relationship). This property is fundamental for bone fracture prevention and for the development of bio-inspired micro-nanomaterials. Osteonal lamellar bone thus displays specific anisotropic values by modifying toughness and strength through the modulation of sublamellar fibril orientation with respect to the loading direction. Both collagen and mineral orientation within the collagen fibril are thus important factors affecting the mechanics of this multi-faceted and hierarchically structured natural composite.