Mechanics – Two way painting algorithm

In our model the structure of trabecular bone is mapped onto a quadratic (2D) or simple cubic (3D) lattice, respectively (lattice constant a and S_x, S_y and S_z, respectively, lattice sites in the specified direction). Each lattice site can be occupied, i.e. filled with bone matrix, or empty, i.e. filled with marrow. Since a vertebra consists of a shell of compact bone enclosing the core of trabecular bone, in the model fixed boundary conditions were chosen

7.1. Mechanics – Two way painting algorithm 67

such that

σ_xyz ≡1 ∀(x, y, z)∈∂V (7.1) with σ_xyz the spin variable at location (x, y, z) that takes the value 1 if the element is occupied and 0 if it is empty. ∂V denotes the boundary of the sample volume V. The architecture of bone is then characterized by the value of the spin variable σ_xyz ∈ {0,1} for each lattice site.

Given a special structure and an external load the local stresses and strains have to be calculated for each bone element. It is assumed that a constant force acts on the vertebra in z-direction. In a first step the load bearing elements of the structure have to be determined in each direction.

It is assumed that the force “flows” like wet paint through the occupied elements Parkinson et al. [1997]. Starting in the plane z = 0 all occupied elements are labelled connected. Then the plane z = 1 is chosen. In a first step all occupied elements directly above a connected one are labelled connected. Then occupied elements to the left and to the right of connected elements are also labelled connected. This procedure is repeated for each layer until the last z-layer is reached. After that the procedure is repeated, now bottom-up, and exchanging the properties connected with active and occupied with connected. After the active elements have been determined in z-direction they are also determined inx- and (in the case of 3D)y-direction.

The result of these calculations is a loaded skeleton, which consists of all elements that are active in a special direction and therefore also experience load and contribute to the load bearing in that direction. Consequently all other elements (empty and non-active) are unloaded. We have

active⇒connected⇒occupied (7.2) That basic version of the painting algorithm leads to an angle of the force flow of 45^◦ (tanα = 1), which is an arbitrary input into the mechanical description. First simulation results suggested that this value is too large to give realistic bone architectures, so the angle of force flow was restricted to smaller values by posing one additional condition: an element to the left or to the right of one connected (active) element was only then labelled connected (active), when the element two layers above (below) the main element was first: connected (active) and second: has a connected (active) element directly above (below). This additional condition leads to an angle of the force flow of ≈18.43^◦ (tanα= 1/3).

Once the active elements have been determined in each direction the local stresses and strains of each of the active elements have to be calculated. The elastic energy W associated with an external force F_z applied in z-direction

68 7. The Model – Stochastic Remodelling with mechanical feedback

Figure 7.2: Schematic representation of the painting algorithm to determine the active (loaded) sites of a special bone architecture in z-direction. The algorithm starts in the top z-layer (z = 0) by labelling all occupied elements (cyan) as connected (green). In the subsequent layers all occupied elements directly beneath connected elements are also labelled connected, as are occupied neighbours directly to the left and right of connected ones. When the last layer (z =S_z−1) is reached the procedure is repeated now from bottom to top to extinguish dead branches. Now the state variables change from occupied to connected and from connected to active (blue). The active skeleton (the set of all active elements) inz-direction are all elements experiencing load in that special direction.

on the sample can be written as follows W =L_xL_yL_z with L_i =aS_i the length of the sample ini-direction, ∆L_i the elongation in i-direction (∆L_i/L_i = _i the strain in i-direction) and E_i the global elastic (Young’s) modulus ini-direction. The first term describes the energy stored in the sample (the energy density ρ of a material with Young’s modulus E deformed by a strain is given by ρ = ¹₂E²), the second term is the work done by the external forceF_z.

In the presented model it is assumed that F_z is acting on the vertebra solely in z-direction. But due to the inwaisting form of a vertebra such a force will result in effective forces, and consequently also in deformations, in directions perpendicular to the original loading direction (x- andy-direction, respectively). In the following the relation of deformations in the loading directionz and the other directions (xand y) is deduced. It is assumed that the macroscopic form of a vertebra may be described as an arc of a circle with constant length L. A deformation of this structure is realized by bringing together its ends by δz, which causes a maximum shift of the middle of the

7.1. Mechanics – Two way painting algorithm 69

arc by δx (see Figure 7.3). Specifying the arc by its radius R and angle θ

Figure 7.3: An arc of a circle with RadiusRand angle2θis deformed by bringing its ends together byδz, which results in a perpendicular shift ofδx. The deformed structure forms an arc of a circle with radiusR⁰ and angle2θ⁰, such that the length of the arc is constant.

before bending and with R⁰ and θ⁰ after bending it is found that

Rθ = R⁰θ⁰ (7.4)

δx = R⁰(1−cosθ⁰)−R(1−cosθ) (7.5) δz = 2Rsinθ−2R⁰sinθ⁰. (7.6) With θ⁰ =θ+δθ one obtains R⁰ = _θ+δθ^θ R and considering only small defor-mations, i.e. δθ 1, sinθ⁰ ≈ sinθ +δθcosθ and cosθ⁰ ≈ cosθ −δθsinθ, where terms of higher than linear order were discarded. Insertion into the preceding equations yields

δx≈ 1−cosθ−θsinθ

2 (θcosθ−sinθ)δz ≡kδz (7.7) which gives a linear relationship between deformations in z- and x-direction with coupling constant k.

Due to the preceding calculations it is possible to connect the deforma-tions in z-direction with deformations in x- and y-direction, respectively, which gives

∆L_x = ∆L_y =k∆L_z (7.8)

Equation (7.3) can therefore be written W =L_xL_yL_z

70 7. The Model – Stochastic Remodelling with mechanical feedback

Minimising this expression with respect to∆L_z leads to

∆L_z = 1 Describing the loading situation with effective forces in all three directions gives

are determined. It is assumed that the material building up bone is isotropic and homogenous with elastic modulus E. Additionally it is assumed that the force in each layer is shared equally between the existing trabeculae. Let M(z)be the number of trabeculae in a layerz =const.andN^j(z)a² the area of the j-th trabecula (j = 1,· · · , M(z)). The strain in the j-th trabecula at height z is therefore given by

^j_z(z) = F_z^?

M(z)a²N^j(z)E. (7.14)

Now the strain is averaged first for the plane z = const. and then for the whole vertebra

Comparing equations (7.13) and (7.15) yields E_z = Sz