A fully implicit and thermodynamically consistent finite element framework for bone remodeling simulations

This work addresses the thermodynamically consistent formulation of bone remodeling as a fully implicit finite element material model. To this end, bone remodeling is described in the framework of thermodynamics for open systems resulting in a thermodynamically consistent constitutive law. In close analogy to elastoplastic material modeling, the constitutive equations are implicitly integrated in time and incorporated into a finite element weak form. A consistent linearization scheme is provided for the subsequent incremental non-linear boundary value problem, resulting in a computationally efficient description of bone remodeling. The presented model is suitable for implementation in any standard finite element framework with quadratic or higher-order element types. Two numerical examples in three dimensions are shown as proof of the efficiency of the proposed method.


Introduction
Bone adapts to the loads under which it is placed. This phrase was stated in 1892 by Julius Wolff and is today referred to as Wolff 's law [41]. In the following, many more researchers contributed to the theory and experimental validation of bone remodeling, among them [16,34], for examples. Based on these insights, connecting mechanical influences and biochemical reactions, many improvements in orthopedic treatment could be introduced into clinical practice. Cowin and Hegedus [11] should be stated here as the first work to provide a closed mathematical description of bone remodeling. This was the starting point for bone remodeling simulations using finite elements, which has continued to be an active field of research in computational biomechanics throughout the last four decades. Despite promising results of even very early models [2,9,40], numerical instabilities, like the occurrence of checkerboard-patterns, were a recurring problem. Without stabilisation, those formulations B Maximilian Bittens maximilian.bittens@ibnm.uni-hannover.de 1 Institut für Baumechanik und Numerische Mechanik, Leibniz Universität Hannover, Appelstraße 9A, 30167 Hannover, Germany were strongly mesh-dependent; refinement results in finer structures with subsequently smaller areas either adopting a prescribed minimum or maximum value for the bone mineral density. Many different approaches have been tried to achieve stability. In [22], for example, a node-based method was introduced, which successfully suppressed the checkerboard modes for linear elements with the density held constant per volume by averaging, which is related to the superconvergent patch recovery method [44]. Finally, Harrigan and Hamilton [19] showed under which conditions bone remodeling with E--relations (cf. [10]) of the form have a stable and unique solution. By this time, three principle approaches for the remodeling stimulus were available in the literature: (1) the stress approach, (2) the fatigue damage approach, and (3) the strain energy density approach. Generalizing these ideas, [8] (cf. [38]) defined a daily remodeling stimulus with N different daily loading cases i, each with n i repetitions and associated stimulus Ψ i , and constants K and m. Choosing an appropriate stimulus Ψ i the approaches (1), (2), and (3) can be recovered. Furthermore, it was shown that if Ψ is uniform in the bone, the three approaches lead to the same basic result. For more information, the reader is referred to the excellent review article [38]. From the more recent past, two different approaches should be mentioned here as examples: [24] and [17]. The former models bone adaption with a strain energy density stimulus built upon a profound theory for thermodynamically open systems, cf. [14,25], and the latter describes bone remodeling in analogy to damage mechanics. It is remarkable that even in the most recent contributions, the most prominent approach for time integration of the constitutive equations is the so-called "staggered"-approach (cf. [21,26,31,35,37]). This approach can be roughly defined by the following three steps: (1) solve a purely linear finite element model of the bone to obtain the mechanical stimulus at the nodes, (2) update the bone mineral densities at each node according to the discretized constitutive function and (3) update Young's modulus of each integration point according to the E--relation. This procedure is repeated until convergence is achieved. The advantage of this approach is that no modifications of any finite element routines are necessary. The major disadvantage is the poor computational efficiency since explicit time-integration of the constitutive equations has to be performed.
It should be noted that there exists a variety of different approaches to achieve stability in computational bone remodeling, for example, [7,15], or [13]. However, all of the latter require some smoothing, averaging, or projection, making them difficult or impossible to implement as a standard finite element material model at the integration point level.
In this contribution, the integration point approach presented in [24] is adopted, but all steps are derived in close analogy to [12], which results in a thermodynamically consistent description of bone remodeling as an infinitesimal strain finite element material model with superior computational efficiency due to the fully implicit formulation. As a novel result, the inner local Newton iteration present in the integration point-based approach in [24] can be omitted.

Constitutive theory
The following briefly summarizes some essential relations between thermodynamically open systems and constitutive modeling. For a complete review of the theory of thermodynamically open systems, the reader is referred to [25] or [23]. In addition, it should be noted here that there exist a set of fundamental principles of material theory, which has to be fulfilled by any constitutive model, for example, the principle of determinism or the principle of objectivity. The reader is referred to [12] or [33] for a complete description of constitutive modeling.

Balance of mass
The local version of the balance of mass for open systems can be stated as the equality of the rate of chance of the spatial mass density and a mass source R, which is left to be defined in Sect. 4. Many other bone tissue models with a balance of mass of type (3) can be found in the literature by, amongst others, [3,19,40]. The incorporation of a mass flux into the latter equation, as shown in [25], is omitted here. In doing so, the resulting set of governing equations would require a numerical discretization scheme, the solution of which would be much more costly. To the authors, this additional expenditure does not seem to be justified since for bone remodeling, both approaches lead to the same basic results, as shown in [24].

Dissipation inequality
Respecting the non-constant mass in an open system results in an additional entropy source S, as shown in [29,36], or [14], for example. Against this background, [25] provided a free-energy density-based version of the Clausius-Duhem inequality for open systems of the form where ε is the linearized strain tensor, the mass density, ψ the specific Helmholtz free energy function, S the entropy, θ the absolute temperature, and Q the heat flux. The Clausius-Duhem inequality can be decomposed into a local term d loc , typically referred to as Clausius-Planck inequality, and a conductive term d cond , typically referred to as Fourier inequality. Both terms are required to hold separately: Here, the assumption is made that all processes are modeled as isothermal processes. Due to that assumption the Clausius-Duhem inequality reduces itself to the local part d loc , since d con ≥ 0 holds if ∇ X ln θ = 0 .

A thermodynamical consistent constitutive law for bone remodeling
Following [12], it will be assumed that the set of state variables determines the thermodynamical state for any time t at a point X ∈ B 0 , where ε is the linearized strain tensor and is reinterpreted as the bone mass density. Consequently, the specific Helmholtz free energy ψ = ψ( , ε) is dependent on the state variables. Using the relation (3), the material time derivative of the free-energy follows as Inserting the above into the Clausius-Planck inequality (5) 1 yields from which the constitutive equations are implied. By that procedure the fulfillment of the Clausius-Duhem inequality is guaranteed a priori. Concluding the above, the thermodynamically consistent constitutive law for isothermal bone remodeling can be stated as Since the constitutive model described above is only dependent on the history of the linearized strain and bone mineral density, it is possible to define the constitutive initial value problem: presuming the history of linearized strain ε(t), t ∈ [t 0 , T ], and the initial value of the bone mineral density (t 0 ) are known, find the history of σ (t) and (t) such that the constitutive equations hold for every t ∈ [t (0) , T ] (c.f. [12]).

Density-weighted generalized Hooke's law
For isothermal processes the strain energy density function Ψ ( , ε) = ψ( , ε) (12) is defined by the product of the bone mineral density and the specific Helmholtz free energy ψ (see e.g. [42]). In bone remodeling, a quite common choice for the specific Helmholtz free energy is based on a classical linear-elastic-type free energy function ψ LE weighted by a relative density ( / 0 ) n , where the exponent n is typically varied between 1 ≤ n ≤ 3.5, and λ and μ being Lamé constants, as shown in [10,18] or [25], for examples. This provides a redefinition of the Cauchy stress tensor as the density-weighted Cauchy stress tensor The density-weighted material tensor can then be derived as Finally, analogous to the generalized Hooke's law for continuous media, its density-weighted counterpart can be stated as

Relation of bone density to mechanical properties
In numerous works, a relation between the bone mineral density and Young's modulus of bone, with the general form of equation (1), has been established (see, e.g. [10,18,30,32]).
Here, E 0 , 0 , and n are left to be identified by experimental investigations and physical reasoning, an ongoing issue in the scientific community. In [27] following basic material properties for the E--relation have been proposed: It has been shown that this model fits experimental observations sufficiently and is therefore used here.

Finite element modeling
This section provides a finite element modeling approach for bone remodeling. Herein, great care was taken on the linearization of the weak form and the constitutive relations and their incorporation into a Newton-Raphson scheme.

Weak form
A strong form of the governing equations describing bone remodeling in the frame of small deformations and with respect to an initial configuration of a continuum body B 0 can be stated as Furthermore, the existence and uniqueness of a solution (cf. [39]) to the strong form (18) is only guaranteed with a suitable set of boundary conditions prescribed on ∂B 0 , namely Dirichlet boundary conditions on Γ u ⊆ ∂B 0 and Neumann boundary conditions on Γ σ = ∂B 0 \Γ u : Equation (19) is equivalent to demanding u to be a kinematically admissible displacement field and σ to be a statically admissible stress field (cf. [39]).
Applying the principle of virtual work and some mathematical manipulations results in a standard weak form for small strain elasticity

Material non-linearities
Presuming a constant material tensor C, the weak form (20) is a linear boundary value problem that in a finite element framework ultimately results in a linear system of algebraic equations. Obviously, the latter does not hold here since in bone remodeling C = C( ) is a density-weighted material tensor (see Sect. 2.4). The evolution of bone mineral density is described by the constitutive initial value problem given in equation (11), a constraint that has to be fulfilled in addition to the weak form. Therefore, the problem becomes a nonlinear initial boundary value problem: with the prescribed history of Dirichlet boundary conditions the prescribed history of Neumann boundary conditions and the initial internal bone mineral density field In general, the constitutive model is path-dependent, and for problem (21), a closed-form solution is not available. Choosing a backward Euler numerical integration scheme for the constitutive initial value problem (11) results in the definition of an incremental constitutive function for the bone mineral density whereR is the algorithmic counterpart of the mass source R and Δt = t (n+1) − t (n) , and an incremental constitutive function for the stress tensor As a next step, the above is reintroduced into the weak form (20), resulting in an incremental boundary value problem. Because within one time-step, the bone mineral density is held constant, we can think of δΠ (n+1) = δΠ(u (n+1) ) as a function of the unknown displacements u (n+1) alone, which makes the constitutive model path-independent within one time-step: The above is a non-linear equation and needs consistent linearization in order to be solved via a Newton-Raphson scheme.

Linearization
Presuming the displacement field u (n) is known and δΠ is sufficiently smooth in t, a expression for the unknown virtual work can be derived at time-step t (n+1) = t (n) + Δt and u (n+1) = u (n) + Δu by truncating a Taylor series expansion: Albeit δΠ only exhibits a variation in u, δΠ was written as a function of and u to indicate the used time-step of each quantity. The same Taylor expansion is applied to the incremental constitutive function, which, by defining Δε = Grad(Δu), can be stated aŝ The variation of the incremental constitutive function can be found by defining the trial strain ε * (n+1) = ε (n) + Δε and applying the chain rule: whereC is the consistent tangent modulus. We can then further define the known principle of virtual work at time-step t (n) : and the increment of the virtual work: with the virtual work of internal forces at time-step t (n) being defined as the virtual work of the external forces at time-step t (n) as and the increments of the internal virtual work and the external virtual work as and respectively, with Δt = t (n+1) − t (n) . Summarizing all of the above, the resultant linearized weak form can be stated as

Discretization in space
In the following, Voigt's notationã for symmetric tensors a is used. Presuming a standard discretization strategy of B 0 with Lagrangian finite elements, a discretized weak form for a generic finite element e can be stated as whereû e is the element displacement vector, H is the element displacement field interpolation matrix and B is the element strain interpolation matrix.û e , and B are chosen such that multiplying results in Voigt notated strains: Note, since Eq. (39) has to hold for arbitrary admissible virtual displacement fields δu, the equivalent system of equations for an element e is written as Subsequently, the assembled system satisfying the mechanical equilibrium equations is written as which can be solved for the unknown incremental displacements Δû, provided valid boundary conditions have been incorporated into the system of equations. For more information on finite element modeling and spatial discretization, the reader is referred to [43] or [6].

Non-linear solution: the Newton-Raphson scheme
In order to solve the non-linear equation (27), the iterative Newton-Raphson method is employed. The Newton-Raphson iteration counter will be denoted by a superscript (k) in parentheses, with the first iteration starting at k = 1, while global time increments will be denoted by a subscript (t) in parentheses.
For the first global time-step, the global displacement vector and the bone mineral density have to be initialized: where 0 is an initial bone mineral density field that can be chosen arbitrarily as long as it does not contradict any of the requirements made above. Note that the displacement field is initialized at the nodal positions while the bone mineral density is initialized at the integration points. Presuming, the solutionû (n) is known, the updated solution vector is computed from the increments Δû (k) of a converged Newton-Raphson procedure with (m) iterations. An incremental version of the weak form (39) Δû e (k) In the above, the algorithmic consistent material tangent can be derived as where the iterative trial strains ε tr,(k) are defined at the integration points, and the stresses are given bỹ Note that the increment of the bone mineral density still only depends on the last converged bone mass density (n) and the trial strainsε tr,(k) (n+1) . Assembling the element weak form (45) results in the system of equations which is solved for the unknown increment Δû (k) . The Newton-Raphson iterations are repeated until, for some timestep (m), the procedure is said to be converged, if the following condition is fulfilled: where tol is a user-defined parameter. Once convergence is achieved, the bone mineral density is updated: Algorithm 1 Material subroutine 1: procedure response(ε (k) (n+1) , (n) ) 2:

Strain energy density-driven bone remodeling
In Sect. 2, a thermodynamic consistent constitutive law describing bone remodeling has been proposed. The balance of mass has been defined in Eq. (3), while the mass source R was left to be defined. According to the principle of thermodynamic determinism R = R( , ε) has to be a function of the state variables { , ε}. In this work, a strain energy density-driven bone remodeling approach is adopted. In [3] a strain energy density approach of the form has been introduced, where is the density-weighted strain energy density for a linear elastic material restricted to small deformations [see Eq. (13)], Ψ ref is a physiological target value that should be adopted by the density-weighted strain energy density, and c is an additional parameter with the unit time divided by area, which governs the speed of the bone remodeling process. In [19] a extension of this approach has been suggested by the introduction of an additional factor ( / 0 ) −m . By setting m = 0 the approach of [3] is recovered, while it has been shown that by choosing m > n, uniqueness and stability of the solution is guaranteed [19]. The necessity of the stability criterion m > n for strain energy density-driven bone remodeling, as introduced in Eq. (55), can be demonstrated by the following example: For an arbitrary point in the continuum we choose 0 = 1 g /cm 3 , n = 2, Ψ LE = ψ LE = 0.1 g /cm 3 , Ψ ref = 0.05 g /cm 3 , and Δt · c = 1 s 2 /m. Since Ψ ref is to be adapted by Ψ LE , it is clear that the actual bone mass density must increase. In Fig. 1, the change of density Δ is plotted for m = 0 and m = 3 over all possible actual mass densities ∈ [0.1, 2]. For m = 0, it becomes visible that, contrary to the physiologically desired behavior, bone mass density decreases for small actual values of bone mass density. As the actual density gets larger, the change in bone mass density increases quadratically. In contrast, for m = 3, the largest change in bone mass density occurs with small actual values of bone mass density and decreases as the actual value increases, which can be seen as a saturation effect. While the former example promotes checkerboard patterns due to the unphysiological behavior, the latter suppresses these patterns successfully.

Material subroutine
During the integration of the local stiffness matrix and the last converged state of the bone mineral density (n) = (X I , t (n) ) as input arguments. By recalling that within the Newton-Raphson iteration of a time-step [t (n) , t (n+1) ], only the trial strains are allowed to be altered (cf. [12]), it gets clear that for obtaining the iterative change in bone mineral density only an evaluation of the function R is necessary since the quantities (n) andε (k) (n+1) are known, and no internal Newton-scheme is necessary. The material subroutine is briefly summarized in Algorithm 1.

Principle of static-equivialent forces and related biomechanical-equilibrated bone-mineral density distribution
At this point, it shall be noted that bone remodeling is a long-term process that takes place over a period of years. This justifies the omission of dynamic forces if the quantity of interest is a biomechanical-equilibrated bone-mineral density distribution [11,20,28]. Now consider a quasi-static example for a linearizedweak form of type (38) for a non-linear but time-independent material: if the surface loads are held constant between two time-steps t (n) and t (n+1) , no Newton-Raphson iteration would take place since the external forces f ext (n+1) are already in balance with the internal force f int (n) . Contrary, in bone remodeling as described here, two (pseudo)-time constants were introduced: (1) Δt for the implicit Euler time integration and (2) c as a constant describing the process speed in equation (53). That results in a possibly out-of-balance right-hand side, although forces are held constant between two time-steps. By that, it is possible to define a biomechanical-equilibrated bone-mineral density distribution: assume t = t(X) and = (X) are given surface loads and a given bone mineral density distribution, It is noted here that, albeit being unconditionally stable for many problems, the time discretization constant Δt in the backward Euler method cannot be chosen arbitrarily large in order for the Newton-Raphson procedure to converge.
What remains is a meaningful definition of t: in [28], surface loads were computed by solving the inverse optimization problem min t 1 2 where ref (X I ) has been obtained by projecting 3D-CT data from a human femur to an associated finite element mesh.

Numerical examples
In this section, numerical examples demonstrate the correct implementation of the bone remodeling algorithm. The entire project was implemented independently in the language [4]. Neither commercial nor open-source finite element software has been used. MUMPS [1] has been used as a solver  For both examples following assumptions are made: to define the bone material, the parameters introduced in (17) are used. The shear modulus is set to ν = 0.3, and a combined parameter Δt · c = 50 s 2 /m is used for controlling the speed of the adaption process. The force is constant and the maximum number of global time-steps is (n max ) = 120. The algorithm is said to be converged if, for some timestep 1 < k ≤ n: | (k) (X I ) − (k−1) (X I )| < tol for all integration points X I . All results are transferred to X-ray images according to [5] since the internal bone mineral distribution is easier to infer from a reduced representation. As a result no smoothing of the results from post-processing has taken place. Note that no color scale is given as this is an arbitrary result of the calibration of the X-ray simulation attenuation law. For a rigorous calibration, one would need several X-ray images of different specimens produced with a standardized X-ray device. Since this was not done, the results are only a demonstration of the performance of the presented method.
All simulation start with a homogenous bone mineral distribution with (0) (X I ) = 1 g /cm 3 for all integration points X I . It shall be mentioned here that the simulation converges to the same results, even if the bone mineral density is modeled as an appropriate smooth Gaussian random field. For this case, the number of external load steps necessary increases significantly. As a remark, it can be stated that in the convergent case, the increment of the bone mineral density field ||Δ || = || (n) − (n−1) || decreases monotonically over the load steps for a constant load vector t (n+1) = t (n) = t.

Model 1: thin plate
As a first example, a thin three-dimensional plate with dimensions of 20 cm×20 cm×1 cm is clamped on its left-hand side, and a shear force q y (y) = 1 kN/mm at 9 cm ≤ y ≤ 11 cm is applied at the right-hand side as displayed in Fig. 2. The reference strain energy density is set to The exponent m is varied in order to study the influence of this parameter. The results of X-ray simulations of the thin plate following the bone remodeling process are shown in Fig. 3 for linear shape functions and in Fig. 5 for quadratic shape functions. Finite elements with linear shape functions. The thin plate model is discretized by linear four-node tetrahedral elements, resulting in a mesh with 709 elements and 816 degrees of freedom by choosing a coarse discretization scheme or, for a finer discretization scheme, in a mesh with 3081 elements and 3267 degrees of freedom.
For the exponent m = 1, it can be seen that strong checkerboard patterns occur for both discretizations, the coarse one shown in Fig. 3a and the finer one displayed in Fig. 3b. Therefore, convergence of the global algorithm is not achieved, while quadratic convergence of the Newton-Raphson method is preserved in almost every load step, as shown in Fig. 4. It should be noted that the quadratic convergence rate is missed when too large load steps are applied, which is probably the case at the beginning of each simulation or when the residual becomes too small due to the numerical tangent modulus used. As the exponent m increases, the bone mineral density distribution gets smoother, as seen by comparing Fig. 3c, e or Fig. 3d,   simulations, without further treatment of these effects, with the methods described here. Quadratic finite elements The thin plate model is discretized by quadratic ten-node tetrahedral elements, resulting in a mesh with 709 elements and 4566 degrees of freedom by choosing a coarse discretization scheme or, if a finer discretization scheme is applied, in a mesh with 3081 elements and 19,032 degrees of freedom. For simulations with m = 1 (not depicted), convergence of the global algorithm was not achieved, and unphysical patterns could be seen. For m = 2, the same unphysical patterns occur but are less pronounced, as seen in Fig. 5a or b. The exponent being m = 3 is the first case where no unphysical patterns are visible, and the global algorithm converges in 13 steps for the coarse mesh and 16 steps for the finer mesh. In Fig. 5c, disturbances can be seen near the support and force application, but these phenomena disappear if the mesh is refined, as seen in Fig. 5b. It can be observed that while m increases, global convergence is enhanced in general. The Newton-Raphson scheme preserved quadratic convergence as discussed for linear elements in the previous example.

Example 2: human femur
A human femur serves as a second example. Boundary and loading conditions were straightforwardly adopted from [27], where loading conditions for the quasi-static case were calculated by solving the inverse problem stated in (59). A total net force of 382 N was applied to the femur. The distribution of the applied forces can be seen in Fig. 6. This static load-equivalent leads to a biomechanical-equilibrated bone-mineral density distribution (cf. Sect. 5.1), which approximately corresponds to the psychological distribution. For more information on boundary and loading conditions for the quasi-static case and the solution of the inverse problem, the reader is referred to [27,28]. Since linear finite elements did not perform sufficiently in the last example, the model is meshed solely by ten-node quadratic tetrahedral elements, resulting in 21,451 elements and 102,048 degrees of freedom. The exponent m is set to m = 4 since this setting resulted in the fastest convergence rate in the last example. At the same time, the reference strain energy density is altered to study its influence. The results are transferred to X-ray images and depicted in Fig. 7. In Fig. 7a, it can be seen that Ψ ref = 0.00001 N /mm 2 is too small as a reference strain energy density and leads to a bone mineral density distribution that is nearly uniform and doesn't develop visible zones of compact and cancellous bone. Increasing Ψ ref leads to bone mineral density distribution which can be considered more realistic, as seen in Fig. 7b or, with even more prominent developed zones with compact and cancellous bone, in Fig. 7c. In Fig. 7d, it can be concluded that Ψ ref = 0.001 N /mm 2 is too large as a reference strain energy with regard to the force applied since, especially in the region of the femoral head, the bone mineral density seems to be underdeveloped. By that, it can be concluded that 0.0005 ≤ Ψ ref < 0.001 can be considered an optimal reference strain energy with respect to the given model, boundary, and loading conditions. It should be noted that these are only preliminary evaluations. For a proper evaluation, a comparison with medical studies would be necessary, which was not done here. Finally, it can be stated that for this example, the Newton-Raphson algorithm preserved quadratic convergence as described for the previous models, and the global algorithm has converged in less than 15 time-steps. The computation time can be specified with less than 10 min on a desktop PC with Intel Core i7 8700k with 32 GB DDR 4 RAM.

Conclusion
In this work, a fully implicit finite element material model for bone remodeling was provided. The theory of open systems was used to derive a thermodynamically consistent description of bone remodeling in analogy to [24]. In delimitation thereof, incorporating the constitutive relations into the finite Bone remodeling followed by X-ray simulation of a human femur with different reference strain energy densities Ψ ref element framework was done analogously to material modeling, as described in [12]. This ultimately results in a novel description of bone remodeling as a fully implicit and consistent finite element material model. Due to the fully implicit description, it is possible to implement bone remodeling very efficiently into any standard finite element framework with shape functions of quadratic or higher order. The additional computational costs caused by the quadratic shape functions can be easily outweighed by the ability of quadratic finite elements to model rounded objects by many fewer elements. The functionality of the described method was demonstrated with two numerical examples. The unsuitability of linear shape functions, with the methods used in this work, can be explained by the resulting element-wise constant strains. This results in element-wise constant strain energy densities, which are then adopted by the bone mass density. It is noted that these and other unphysical patterns, as depicted in Figs. 3e, f and 5c, d, would not be readily visible from a standard finite element postprocessing. From that, it can be concluded that the results shown in Fig. 7 are smooth on the integration point level without any need for further numerical treatment.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.