Torsion - an underestimated form shaping entity in bone adaptation?

OBJECTIVES: There is ample agreement that the specific shape of a bone is related to the loads it has to carry. It is also believed that bones mechano-adapt in order to 'find' this shape. The open question is which signals constitute the determinants of this adapation. Recent in vivo experiments show that torsion is a significant load component in human tibia, and a computational study of the mechanostat has indicated that torsion could play a role in the shaping of tubular long bones.
METHODS: An earlier computational approach is further progressed to systematically study the relative importance of axial compression, lateral bending and axial torsion.
RESULTS: Results demonstrate that shape-driving potential towards tubular shapes is greatest for torsion, followed by bending and least for axial compression. Multiple linear regression analysis confirmed the dominant role of torsion, in particular for the 2nd moment of intertia. The obtained results were largely unaffected by starting conditions, e.g. either from a grid or through reshaping under disuse.
CONCLUSIONS: Strong support has been found for the hypothesis torsion could be more important than suggested in previous studies as a component of the mechanical environment of bones. This will apply to the shafts of long bones, and also to the femoral neck.

Introduction

Bone is a living tissue that is continuously remodelled[1-3], normally without obvious modifications of macroscopic structure and shape. This changes when significant alterations in the loading spectrum occur, be it because of disuse due to reduced gravity in space or long bed rest[4-6], during the recovery after bed rest of space-flight[7,8], or in the course of sport-specific increased loading[9-11]. In such cases we observe significant modifications of bone geometry showing an increased bone mass with increasing load and reduced bone mass with decreasing load. To describe the underlying regulatory mechanisms in a formal way, Frost proposed the ‘mechanostat’ model[12,13]. Such a formal description of bone mechano-adapation has been the basis for clinical differentiation of primary and secondary bone disorders[14,15].

A plethora of attempts exists to elucidate the rules of bone mechano-adaptation, starting with the qualitative theory known as Wolff’s law[16], over Frost’s semi-quantitative three-way-rule[12,13] up to more elaborate in silico studies (for survey see[17]). Of the latter, the model proposed by Huiskes and co-workers[18,19] is particularly attractive because it differentiates between dynamic and static loading[20], and also because it is biologically-motivated, namely having the osteocytes, living cells, embedded inside the matrix of the compact bones, working as mechano-sensors[21].

Recently, we have applied that approach to perform a ‘computer-assisted experiment of thought’ into the shape forming processes of shaft-like geometries[17]. Results of that study yielded, as expected, the sensitivity to model parameters, but additionally and quite surprisingly that the Huiskes mechanostat can explain flexure neutralization, the adaptive straightening of flexed children’s bone after ill-healed fractures. This latter result was in stark contrast to the expectations in literature[1]. Another relevant finding was that torsion might be more important than suggested by previous studies[22-26]. Moreover, recent work has found torsion to be a significant component in the load spectrum of long bones in vivo[27].

Therefore, the current study has been designed to further clarify this latter aspect, and to disentangle the relative influence of the different load types on the geometry-shaping process. Accordingly, the loads imposed by axial compression, bending and axial torsion, individually or in combination, were systematically varied, in order to study the induced alterations on our standardized starting geometries.

Methods

Transformation cycle

The mechanostat algorithm used in this study has been described before[16].

In brief, each transformation cycle (section B of Figure 1) starts from a given model geometry to which a random load pattern comparable to those acting on the diaphysis of long bones is applied. Next, the resulting stress distribution within the structure is calculated by means of Finite Element Analysis (FEA) using commercial software (ANSYS version 12.1; ANSYS Germany GmbH, Darmstadt). Strain energy density rate in the volume elements of the model serves as a control parameter and governs the transformation algorithm, which in turn modifies the model geometry in an iterative way. The in silico design space was a cube of 4.5 mm edge length, consisting of 421,875 cubic voxels with 60 µm edge length each. “Bone” and “no-bone” voxels were implemented by the setting of the related material properties, set to 15 GPa respectively 200 MPa for bone or no-bone voxels. The Poisson number was set to 0.36.

Figure 1

Part B shows a sketch of the transformation cycle; further explanation found in the text. Part A shows the geometries that have been chosen as starting geometries; the left panel of the image shows the cylindrical lattice together with the related cross-section silhouettes (horizontal mid-shaft and vertical) where the block 1 simulation of CL type started; the right panel shows the straight closed tube that itself is the product of our simulations together with cross-sectional silhouettes, used as starting geometry for the FT type simulation runs. Part C represents a sketch of the loading pattern; the black arrow represents a force on the rigid upper plane with a constant axial component in z-direction combined with a randomly distributed lateral component in the x-y-plane; the magenta double head arrow represents a torque on the plane around its centre of area. All three sketches appear in similar form in [17].

Starting geometries and their generation

Two different brick models were used as starting geometries, defined by the setup for bone properties in the corresponding voxels. In z-direction, the geometry extended over the whole width of the design space, in x and y directions the diameters were smaller than the design space leaving space for the evolution of the geometry. Two starting geometries were used in order to investigate a) the building of a structure from a neutral lattice-like precursor structure and b) the maintenance (or degradation) of a given bone–like geometry. A depiction of the starting geometries can be seen in section A of [Figure 1]:

a cylindrical shaped cubic Bravais lattice model with the axial length of the full design space size of 4.5 mm and a diameter of 3.38 mm. The unit cell was a cube with an edge length of 8 voxels (=0.48 mm) with beams of 2x2 voxels at all edges building a cubic frame. Combining these unit cells, a cubic lattice of quadratic beams with a diameter of 4 voxels (=0.24 mm) was generated. For the creation of the CL we used a template lattice generated in Matlab (Mathworks Inc., Natick, Mass. USA version R2011a) and matched that lattice into the mesh of our FE model using ANSYS standard functionality. Irregular surfaces of the lattice rods were due to this matching process. The irregularity was assumed to not influence the simulation results, because the same two starting lattices were used for all simulation runs. The CL was the main starting geometry for the investigation.

In order to investigate the reversibility of effects, a total ‘bone disuse’ was modelled using an – up to then – fully maintained tubular structure. Note that in the context of this study “full tube” signifies a hollow cylinder with “fully” closed walls. This starting geometry derived from a previous simulation run with 100 N axial compression, 40 N bending load and 180 Nmm torsion, resulting in a regular tube.

The loading pattern

The principal loading pattern is shown in section C of [Figure 1]. We applied axial compression forces of up to 300 N, lateral bending forces of up to 40 N and axial torques of up to 180 Nmm in different magnitude levels and different combinations of those levels. Note that the bottom of the structure was always fixed. It should be noted as well that due to technical limitations, dimensions had to be downscaled in comparison to human scales. For the justification of the downscaling and discussion of its potential effects we refer to our previous study[17]. For the selection of the loading magnitudes and the scaling problem when comparing our mouse size model with human geometries see appendix A.

Transformation law and its implementation

The transformation law has been adopted from Huiskes et al.[18,19] and simulates the modification of bone density m (value between 0 and 1) in terms of the bone-generating activities of the osteoblasts and the bone-resorbing activities of osteoclasts. Following Huiskes et al., strain energy density (SED) rate was taken as the leading mechanical signal that activates the osteocytes. We refer to our previous publication[17] for a detailed description of the algorithm.

The transformation cycle was implemented in Matlab. Strain energy densities resulting from the simulated load applications were calculated using FEA based on the commercial package ANSYS, using a Preconditioned Conjugate Gradient (PCG) solver. Calculations were performed on a high performance computing cluster (IBM x3550 and X iDataPlex computers). For each run, between 200 and 400 iterations were calculated to make sure that the forming process did converge. Each calculation step (a single loop of the transformation cycle) took up to 60 minutes of computation time, depending on the number of remaining bone voxels. Considering an average step-duration of 30 minutes and a total number of 56 runs with 200 steps each we arrive at an estimate of 240 days of computation time in total.

Simulations overview

We performed two blocks of simulation:

In the first block we implemented any combination of three magnitudes of load w.r.t the three load types axial compression, bending and axial torsion. In these calculations we started with the CL geometry.

1) Axial compression force in z-direction was applied as a constant value with the following magnitudes:

a) No compression force, b) 50 N, c) 100 N, d) 300 N.

2) Bending force was applied on the top surface in an x-y-direction randomly selected over 360° with following magnitudes:

a) No bending force, b) 20 N, c) 40 N.

The bending arm length being 4.5 mm, that settings correspond to bending moments of 0, 90 and 180 Nmm. In the present study we do not consider directed bending. Note that the fixed bottom generates a cantilever type of bending that has to be taken into account when comparing with real biomechanical situations.

3) Axial torque was applied around the central z-axis of the top surface as a random magnitude between:

a) No torque, b) 90 Nmm ccw and 90 Nmm cw, c) 180 Nmm ccw and 180 Nmm cw.

Thus, a total of 36 computational runs have been performed for that block.

Note that on a human scale those 300 N axial compression, 180 Nmm bending moment and axial torque of 180 Nmm as well would correspond to 10.4 kN compression, 39 Nm bending moment (bending load of 1440 N on a 27 cm arm) and 39 Nm axial torsional load respectively. For a deeper explanation refer to appendix A. In appendix A we are also showing that the maximum stress generated in the bone (for the ideal geometrie used there) for our maximum load cases are the same for all the load modi. For our simulation we extrapolate a maximum stress level of 30 Mpa for all modes of load.

In the second block we started with the FT geometry. In total 20 additional runs have been performed for all bending and torsion cases, with the axial compression restricted to zero or a 100 N load.

An overview of all runs is given in [Table 1].

Table 1

Overview of the computational runs performed in this study. ID signifies a run identifier that contains the block code (first 2 characters). The numbers in the code represent the different load levels for compression, bending and torque at positions 3, 4 and 6, respectively. Z gives the axial load in N, XY the lateral bending load that was applied with random direction from step to step. Torque is the range of axial torque applied; e.g. ±90 means that within a range from 90 Nmm clock-wise to 90 Nmm counter-clock-wise torques have been applied randomly with uniform distribution. #steps indicates the total number of steps computed for this computational run. Starting geometry indicates one of the two starting geometries as previously described. In all runs we used the following settings for the Huiskes parameters: τ=5⋅10-7 mm5/(Nmol), kthr=0.05, D=100 µm, µ(osteocyte)= 1 nmol·mm·s/(J·day), F=1 Hz, Vr=1.5 mm3, focl=7.1·10-4/(voxel·day), osteocyte density=44⋅103/mm3. See references[16-18] for details of the model.

ID	Z (ax.) [N]	XY (lat.) [N] 360°	Torque [Nmm]	#steps	Starting geometry
CL0000	0	0	0	93	CL
CL0001	0	0	+/-90	194	CL
CL0002	0	0	+/-180	199	CL
CL0100	0	20	0	198	CL
CL0101	0	20	+/-90	300	CL
CL0102	0	20	+/-180	191	CL
CL0200	0	40	0	178	CL
CL0201	0	40	+/-90	191	CL
CL0202	0	40	+/-180	196	CL
CL1000	50	0	0	197	CL
CL1001	50	0	+/-90	196	CL
CL1002	50	0	+/-180	186	CL
CL1100	50	20	0	187	CL
CL1101	50	20	+/-90	191	CL
CL1102	50	20	+/-180	200	CL
CL1200	50	40	0	200	CL
CL1201	50	40	+/-90	200	CL
CL1202	50	40	+/-180	200	CL
CL2000	100	0	0	200	CL
CL2001	100	0	+/-90	200	CL
CL2002	100	0	+/-180	200	CL
CL2100	100	20	0	200	CL
CL2101	100	20	+/-90	200	CL
CL2102	100	20	+/-180	400	CL
CL2200	100	40	0	200	CL
CL2201	100	40	+/-90	400	CL
CL2202	100	40	+/-180	600	CL
CL3000	300	0	0	200	CL
CL3001	300	0	+/-90	200	CL
CL3002	300	0	+/-180	200	CL
CL3100	300	20	0	200	CL
CL3101	300	20	+/-90	200	CL
CL3102	300	20	+/-180	400	CL
CL3200	300	40	0	200	CL
CL3201	300	40	+/-90	400	CL
FT0000	0	0	0	200	FT
FT0001	0	0	+/-90	1200	FT
FT0002	0	0	+/-180	200	FT
FT0100	0	20	0	200	FT
FT0101	0	20	+/-90	200	FT
FT0102	0	20	+/-180	200	FT
FT0200	0	40	0	200	FT
FT0201	0	40	+/-90	200	FT
FT0202	0	40	+/-180	200	FT
FT1000	50	0	0	200	FT
FT1200	50	40	0	200	FT
FT2000	100	0	0	200	FT
FT2001	100	0	+/-90	200	FT
FT2002	100	0	+/-180	200	FT
FT2100	100	20	0	200	FT
FT2101	100	20	+/-90	200	FT
FT2102	100	20	+/-180	400	FT
FT2200	100	40	0	200	FT
FT2201	100	40	+/-90	400	FT
FT2202	100	40	+/-180	200	FT

Criteria for convergence were defined as previously described[17]. In an attempt to classify the geometric and topological properties of the resulting structures, the following parameters have been computed for the x-y-planes (slices) at z-position (axial direction) 38 out of 75:

1) The number of ‘bone’ voxels NoV in a given cross-section.

2) The maximum principle second moments of area MSMoA around the center of that cross-section in the design space.

The focus has been set on the mid-shaft, where the influences of boundary effects of the force injection at the upper and lower surface of the structures are at their minimum.

Finally, Multiple Linear Regression was performed in order to get a quantitative estimate of the relative bone-shaping effects of compression, bending and torsion. To this purpose, compression, bending and torsion data were re-scaled to the range of [0..1], and NoV and MSMoA were z-transformed. The regression analysis thus yielded beta values, i.e. regression coefficients that are normalized to the standard deviation of the dependent variables.

Results

Block 1 results: starting from cylindrical lattice with random bending directions

[Figure 2] provides an overview of the resulting geometries in two different perspectives, one of them (section A) a 3-dimensional representation of the full geometry, the other (section B) a mid-shaft cross-sectional representation. In general, the geometries tend to be more ‘shaft’-like tube structures towards the right (increasing torque). Towards the top (increasing bending) this trend is enhanced. Towards the top (increasing bending) and towards the rear (increasing compression) this trend is increased even further. Of special interest: there is no ‘shaft’-like structure without torsion! A more systematic analysis of these figures reveals the following.

Figure 2

Section A gives an overview about final geometrical results of the 36 block 1 simulation runs, given in form of single 3D bodies that represent the resulting “bone” geometries; The shaded planes have been introduced to help as optical guiding means in the “3D-space” of torsional, axial compression and lateral bending loads. The number of steps in the calculations were approximately 200. The bodies have been generated as 3D-surface models in STL-format directly from the pattern of the simulation space with Matlab. The appropriate scenes of [Figure 4] have been rendered with the open source tool Blender (Blender Foundation, Amsterdam). The representation is organized in a way that three planes with pseudo-3D diagrams are stacked above each other to represent the three lateral bending levels 0, 20 and 40 N. Within the planes, the two remaining dimensions represent axial bending from left to right and the axial compression from front to background. Section B of this figure gives an overview about the transversal cross-sections of all geometries at mid-shaft (50%); representation is analogue to that in [Figure 5]; every cross-section is an image of 75x75 voxels in binary representation with bone in white and no-bone in black. In analogy to [Figure 4], transversal cross-sections at 50% of the bodies are plotted in order to provide additional information on the body topology in that region. The mid shaft area is of particular interest, as it minimizes influences of border effects from the force transfer at the lower and upper plane. As expected, the simulation algorithm produces a no-bone space for an unloaded geometry. Therefore the related position in both figures remains empty.

Figure 4

Final geometries of block 2 calculations; starting geometry was the result of block 1 calculations under high compression, high bending, high torque conditions; only zero or 100N compression force were applied.

Figure 5

Number of bone voxels (NoV) and maximum second moment of area (MSMoA) at 50% for the “remodelling” case; same representation as in [Figure 3].

Compression only (left row of bottom diagram) generates compact rods that increase in diameter with increasing compression force. As demonstrated by the cross-sectional images in section B of [Figure 2], the resulting structure is always solid.

Bending only (shown in the left and front cross-sections of the upper two diagrams) shapes strongly conic forms that are of tripod-like character. With increasing magnitude of bending torque the structure tends to close from top to bottom but with an irregular shape to be seen in the cross-sections of section B in [Figure 2].

Torsion only (bottom row of bottom diagram) generated a tubular truss at 90 Nm and an open-walled cylinder at 180 Nm of torque. Notably, the latter result was still open towards both ends (section A of Figure 2), but actually closed at mid-shaft, showing the typical ring-like cross-section of long bones (section B of Figure 2).

The combination of compression with bending (left outer plane when combining all diagrams) led on all bending levels towards a closing of the geometry untill the emergence of a solid cylinder. Markedly, within this load combination, the effects of compression (towards a solid cylinder) were dominating over the effects of bending (towards a tripod-like geometry).

The combination of compression with torsion (lower plane) quite clearly leads to a central erosion and formation of a central cavity at the highest compression level. Only in its highest level is compression able to reduce the open trusses of the mid-level torsional regime.

The combination of bending with torsion (“front plane” when combining all diagrams) generate an almost perfect tube at the highest simulated bending and compression levels. At lower levels, torsion seems to mitigate the conic and irregular shapes that occurred with bending only. Conversely, bending seems to have potential to close the trusses induced by pure torsional loading.

Thus, none of the three loading types are able to generate a shaft-like geometry by themselves. However, these qualitative results suggest that torsion is a prerequisite for this, albeit in combination with either compression or bending.

“Best” results are achieved when all loads types are present. Notably, there was virtually no tendency by bending, either alone or in combination with compression, to drive the geometry towards a closed cylinder shape.

[Figure 3] takes these qualitative results a step further and provide a quantitative analysis.

Figure 3

Numerical analysis of the cross-sectional images in [Figure 3]; depicted is the commutation of the number of bone voxels (NoV) as well as well as the maximum second moment of area (MSMoA) at the mid-shaft (at 50% of its length); NoV reflects the ‘amount’ of bone available at the investigated cross-section and thus serves as a surrogate for strength in axial compression. The three diagrams A-C show the dependencies of NoV from the perspective of the three different load types axial compression (A), axial torque (B) and lateral bending (C). MSMoA is a measurement for the resistance to bending and torsion. Its dependencies on the load types compression, bending and torsion are shown in plots D-F respectively. The line coding is as follows: Torsion is coded in color: green=0, red=90, blue=180 Nmm; neutral =black; bending is coded in dot size: small=0, middle=20, large=40 N; neutral =cross; compression is coded in line type: dotted=0, small dash=100, large dash=300 N; neutral=full line; for clarity reason we omitted the 50N compression curves for the B, C and E, F diagrams. Note that torsion color coding is ranked as green < red < blue, indicating in the plots A, C, D, F that the number of voxels as well as second moment of area steadily increase with torsion.

It can be seen from [Figure 3], diagrams A-C, that no single load mode dominated the impact on resulting bone strength represented by number of voxels NoV, as evidenced from the scatter of curves in all three diagrams. The maximum value for NoV is generated when all load mode magnitudes are highest. For several load combinations the increase in NoV is quite steep and nearly linear with compression magnitude, for others - those with higher torsion and bending contribution, especially contributions of both - it shows only small increase or even a flat for low magnitudes. Comparing the 0-bending-high-torque curve (blue, small dots) with the high-bending-0-torque curve (green, large dots) in diagram A we find both on a quite similar level, the torsion with light advantages. So we might conclude that torsion fosters bone strength at least to a similar extent as lateral bending. The increase of NoV with torque is moderate and similar to the effect of bending. In maximum we see increases of some 500 voxels looking at absolute values of 400 to 2000 voxels. NoV seems to be very much enhanced by compression when increasing compression levels as can be seen in [Figure 3] in the upper left diagram A. Here we can see jumps of 800 voxels, going from 100 N to 300 N axial compression load.

By the same token, diagrams D-F of [Figure 3] demonstrate a strong dependence of the MSMoA on torsional loading (as evidenced by the narrow overlap of curves in diagram E of Figure 3). The increase of MSMoA ranged from 0.3 mm4 to 3 mm4 for zero-torsion to absolute values between 4.2 and 5.5 mm4 for maximum torque. The dependence of MSMoA on compression and on bending is much less pronounced. Here we see peak increases of 1.5 mm4, yet in certain cases even decreases of 1.5 mm4 e.g. for increasing bending load at high torsion and compression levels. These impressions were quantitatively corroborated by multiple linear regression analysis. For compressive strength (assessed via NoV) we obtained significant correlation (beta coefficient with P<0.001) with compression (1.38±0.16), bending (1.71±0.21), torsion (2.10±0.22), and torsion combined with bending (-1.44±0.33). For the MSMoA significant correlation was found only for torsion (2.10±0.19). All other combinations did not show significant correlations for NoV or MSMoA. The results indicate that compression, bending and torsion had more or less equal contributions to the compressive strength (assessed via NoV), whilst the MSMoA was dominated by torsion alone.

Summarizing, it can be stated that compression has a strong positive effect on the NoV when a certain compression amplitude is exceeded. The highest positive influence on MSMoA has torque. Compared to that, the effect of bending is moderate for both indicators.

Block 2 results: starting from straight closed tube with random bending directions

Instead of starting with a cylindrical lattice the initial geometry was a well-shaped tube (FT). This allowed us to investigate the extent to which the final geometry depends on the starting geometry or, vice-versa, whether the shape-forming could be reversible. The convergence for these calculations took slightly longer than for the CL approach in block 1. In order to verify the results we extended the calculations of two of the runs (FT2201 and FT2102) to 400 steps, but did not find any significant changes. For this principle verification we considered the two cases 0 and 100 N axial compression only.

[Figure 4] depicts the related geometric and topological pattern that we showed earlier in [Figure 2]. Because the processes under observation in this block resemble ‘remodeling’ we will refer to this term in order to distinguish it from the ‘modelling’ in block 1.

In the high torsion regime we find a quite high similarity between ‘modelling’ and ‘remodelling’ results. That can be seen in the outer 3D-view on the geometries as well as from the cross-sections at 50%. In the low torque regime we have minor differences for the 0-bending cases where we have slightly different characteristics in the degree of openness in comparison to the ‘modelling’ experiment.

Larger differences can be found in the 0-torsion regime. Obviously, the different starting conditions lead to multiple rod like structures for the compression only cases and significantly reduces the conic character of the geometries at low and high bending. The structures of the high bending cases are slightly more regular for the ‘remodelling’ case as can be seen from the cross-sections at 50%.

[Figure 5] illustrates the indicators NoV and MSMoA for the ‘remodelling’ experiment series. Whereas we have quite a good compliance for the NoV the MSMoA values do not decrease as much for the 0-torsion case, resulting in a more flat torsion dependence. That correlates with the observation at the geometries themselves that the generation of cone-type geometries in the low torsion regime with bending load elements is reduced and substituted by more truss-type of structures. This results in similar NoV values on the one hand but in significantly higher MSMoA-values on the other hand.

Therefore, it might be the case that some dependence exists on these modified ‘remodelling’ starting conditions, and that this dependence is greatest for low levels of torsion and compression, and is thus mainly related to bending. Obviously it is difficult for the mechanostat to reconstruct a tripod starting from a tube structure. The principal findings of the ‘modelling’ calculations, however, are unaffected.

Discussion

An in silico experiment has been performed in order to find out how the different types of mechanical loading, originating from axial compression, lateral bending and axial torsion, might contribute to build-up and maintain the skeletal structure and its mechanical properties, e.g. that of a tubular long bone. The hypothesis was that torsion might play a considerable more important role in these processes than previously believed. For the calculations a mechanostat model by Huiskes et al. proposed for remodelling of trabecular bone has been adopted. A 3-dimensional loading-space of four compression, three bending and three torsion magnitudes has been analysed. The lateral bending load case has been applied focusing on a uniform distribution with respect to the bending direction. Because we wanted to know whether our results depend on the starting conditions we calculated two different ones: a ‘modelling’ experiment starting from a cylindrical lattice and a reduced ‘remodelling’ program starting from one of the tubular shaped geomatries resulting from the prior ‘modelling’ simulation. Besides the 3D-geometry and the topology of the cross-section at midshaft we surveyed the indicators NoV and MSMoA at mid-shaft as well. We used a model with scaled geometrie with respect to human proportions. In translation our scaled load magnitudes cover a representative range between 0 and 10kN for axial compression, 0 and 29 Nm bending moment and a similar range of 0 to 29 Nm for axial torque. The maximum magnitudes are 30% above known physiological maxima in terms of maximum bone stresses. We have shown that for an ideal tube our maximum loads for all three modes would generate a very similar maximum stress in the material of some 30 MPa, so that the load cases can be assessed as comparable in magnitude. It can be assumed that in principle all kind of combinations of the load types can appear under sports conditions and at reduced magnitudes during normal life. For the bending case it has to be taken into account that cantilever-like bending load on a long bone is the exception because the motion of joints seldom are totally fixed but it can not be completely excluded.

For the ‘modelling’ approach it has been found that there is a broad range of geometries generated ranging from thin compact rods when under axial compression only, over conic tripod-like structures under bending-only conditions and truss-like formations when under low torque loads up to the expected closed tubular shape. The latter is only developing when multiple load types are applied, out of which one has to be axial torsional load. It is not developing when bending and compression alone are applied. Both exotic trends, the conic shaping under bending resulting partially from the cantilever character of our set-up, as well as the multi rod shaping of under compression are “cured” by increasing torsional contributions. In order to reduce effects of the cantilever problem we restricted our further analysis to the central part of the model shaft which we saw as a largely unbiased test bed for our simulation purposes.

The MSMoA indicator significantly demonstrates the highest positive dependence on axial torque. Only the NoV indicator is influenced stronger by axial compression. For the variation of characteristic at lower compression levels from steep to flat we dare the following explanation: increasing contributions of torsional and bending loads, especially in combination, generate sufficient bone strength that compression contributions become effective only behind a certain threshold of magnitude. That indicates the importance of both torsional loads and bending loads. We are quite sure that a systematic underestimation of the compression level might not be the reason.

In the ‘remodelling’ part we found a clear reduction of the structures with respect to closeness of shape and quality indicators. The final results accorded, to a broad range, to the geometries developed by the ‘modelling’ experiment. We found deviations especially in the low torque regime where the conical forms typical for compression plus bending load combinations were substituted by open truss-like forms. This also becomes manifest in a lower dependence of MSMoA on torque while NoV shows the same behavior as in the previous ‘modelling’ case. At the moment we do not have an attestable explanation for that behaviour. It may have to deal with the fact that the thin upper top of the cone structures cannot develop because of the empty space in this area in the tubular starting geometry of the ‘remodelling’ case.

The limitations of the study become obvious if we look on the no-load case. When the simulation is running long enough all bone voxels disappear. Real bone does not disappear when disused for long time and it does not extensively modify its shape. The assumption of bone being resorbed randomly and uniformly distributed needs modification. There might of course be one or more additional control elements that are not load-driven only. Yet, the principle results of this study are not affected by this part of the model. Of more interest might be the question as to “how the choice of the regulating parameter – being the strain energy density as suggested by Huiskes et al. – could influence the structural response not only in magnitude but also in principle?”

Another possibly limitating issue might be addressed as well. In our study we assumed implicitely that maximum loads define the final geometries. We did not check which role the frequency of maximal events might play. Which consequences will it have for instance, if we apply a given magnitude to different directions not in a unique distribution, but with an odd distribution? Will we see an asymmetric final geometry at the end? In computations with explicit directed bending set-ups we saw resulting asymmetries, but in that study the less-loaded directions were never loaded at all. In our torsion set-up we deal as well with a reduced occurrence of maximum load because we admit intermediate torsion magnitudes. The answer to this question might also affect the way of comparison we performed with our indicators above. Can we compare the behaviour of MSMoA for our all-directions bending cases (where one selected direction is loaded less often) with the compression cases where the load is there in every step? If we think further we might expect as well that the asymmetric distribution over the load types in real life - with compression being by far the most often load for a long bone - might bias the result correspondingly. Our answer is that maximum load is the driver – and hypothize that for a future study. Frequency might play a minor role or play a role insofar that a minimum frequency of load cases should occur.

Conclusion

An extensive study on the structural response due to different combinations of compression, bending and torsion acting on bone-shaft-like structures under the regime of the Huiskes mechanostat model has been performed. In the structure forming regime we found a strong support for our hypothesis that the presence of torsional loads is required in building up the well-known tubular structure of long bones. In the structure maintaining regime we found torsion having a strong positive influence in preserving the tubular closed form.