Multiscale investigation of the functional properties of the human femur

1. Introduction

The human skeleton performs some obvious functions: shape and support; attachment of ligaments and muscles; articular leverage in movement; and mechanical protection of vital organs. It also has two vital metabolic functions: haemopoiesis (generation of blood cells), which takes place in red bone marrow, and calcium homeostasis in the blood, which is ensured by controlled dissolution of some mineralized bone matrix during periods of low calcium intake. Primarily, two specialized cell types regulate the maintenance of the mineralized bone tissue: osteoclasts, which destroy the mineralized collagen matrix, and osteoblasts, which produce new collagen that is subsequently mineralized into new bone. The modulation of the replication, activation and apoptosis of these two cell populations, and of the mineralization process, produces the net balance of this metabolic process. In healthy conditions, the amount of bone tissue that is reabsorbed is equal to the amount of newly formed bone tissue, and the total bone mass remains unchanged. Unfortunately, multiple pathological conditions can alter this balance, thus producing a metabolic imbalance of catastrophic consequences.

Three typical examples of diseases, of particular social and clinical relevance, are as follows.

1. Osteoporosis. This produces an alteration of the bone metabolism: the bone mass decreases and the risk of spontaneous bone fractures increases over time.

2. Bone metastases. These are secondary tumours that form in the skeleton as a result of the spread of malignant tumour cells in the blood stream. The formation of bone metastases, besides oncological implications, is a source of tremendous discomfort for these patients. Bone metastases are painful, increase the risk of spontaneous bone fractures, which heal with difficulty and slowly, and can increase the level of calcium in the blood to pathological levels.

3. Renal osteodystrophy. This is a bone pathology characterized by defective mineralization that results from renal disease. It is relatively common in dialysed patients.

These diseases have little in common. The first one is developing as a social pandemic in association with longer life expectancy (NIH 1984, 2000; WHO 1994) and less active lifestyles (Lanyon 1980); the resulting fractures at the hip, the spine and the wrist are very disabling and involve a significant morbidity in the elderly (Lochmüller et al. 1998, 2000; Patel & Murphy 2006). The second one is the result of metastatic tumours; it is rarely a source of death, but is probably the primary cause of pain in oncology patients. The third one is a rare complication (it is frequent only in patients with kidney failure), but it poses extremely complex clinical problems. But the three aforementioned pathologies have one thing in common: they produce a weakening of some portions of the skeleton, which can result in macroscopic or microscopic fractures, even without traumatic events (Rockwood et al. 1991; Coleman 1997; Dennison et al. 2006; Schumock & Sprague 2007).

The primary function of the skeleton is to withstand the biomechanical loads that our daily life produces, without significant damage, internal and external. The propensity of the skeleton to fracture under physiological or paraphysiological but subtraumatic loads (e.g. occasional overloads such as those due to stumbling) is clinically called bone fragility. This term, however, has a certain degree of ambiguity in the bioengineering context. Thus, we prefer to refer to the ability of the skeleton to withstand loading successfully during its daily function as functional competence.

Thus, the study of any means to interfere positively with bone metabolism is vital for the treatment of a number of diseases of great social and clinical impact. Unfortunately, the assessment of efficacy is extremely complex for two reasons. First, each animal model is able to reproduce only certain aspects of the disease and the metabolic processes that the drug alters. A complex combination of various animal models is thus necessary. However, the use of animal models is limited by a non-obvious interaction between factors, although this problem is substantially ignored, even if its importance is well known. That is why bone drugs are always one of the first candidates for in silico assessment projects (Defranoux et al. 2005).

The second reason is that the problem is inherently multiscale and systemic. The drugs operate at the molecule level, their effect is expressed at the cell level and the final clinical effect is produced at the tissue level. The strength of the bone is related to the quantity of bone tissue, its chemical composition and its spatial organization. All three factors are the result of a complex concurrent process at cellular scale. Any drug interferes with a complex cascade of molecular events that modulate the replication, activation and apoptosis of bone cells. Assessing the efficacy of a bone drug requires modelling molecular, cellular and biomechanical processes, each level being entangled with the others. The only way to build effective in silico assessment solutions for bone drugs is to use a multiscale modelling approach (Defranoux et al. 2005).

Although classic reductionist biological research is elucidating more and more complex mechanisms on the genetic and molecular modulation of bone metabolism, so far this knowledge has not translated into a better understanding of skeletal fragility as seen in clinical practice. The approach that Denis Noble and Sydney Brenner called middle-out (Noble 2006) is much more effective: as the event of interest (the fracture) occurs at the organ level, it is wiser to start from this scale, and then drill up, down and across other levels and subsystems.

However, femoral fractures have never been investigated as a multiscale problem. In order to predict the risk that a given patient will suffer a low-energy fracture in the coming years, a multiscale approach is needed. Such risk cannot be addressed at a single scale: the risk of overloading is defined at the body level and the risk of fracture at the organ level. In addition, to predict the evolution of bone strength over time, it is necessary to correlate the molecular cellular processes with the resulting functional effects at the organ level. In this paper, a subject-specific multiscale model is envisioned that goes from the body level down to the cell level. Such a model can be used to generate predictions of the risk of fracture much more accurately than today (most current clinical implementations rely only on bone mass measurements (WHO 1994; Melton et al. 2005)).

This paper exemplifies this approach by describing a model capable of assessing the functional competence of the proximal femur. This example is chosen owing to its clinical relevance (this is one of the most frequent fractures (Dennison et al. 2006)), and for this reason this is the most extensively investigated musculoskeletal region and why a more complete multiscale approach is available. State-of-the-art investigating and modelling techniques are deployed in this example. Phenomenological investigations of fractures have been reported in the past, as well as numerical models (Backman 1957; Cotton et al. 1994; Lang et al. 1997; Keyak et al. 1998; Lochmüller et al. 1998; Keyak 2000; Eckstein et al. 2004; Mayhew et al. 2005). However, femoral fractures have never been investigated as a multiscale problem.

We show how it is now possible to generate, from subject-specific data, an organ-level model able to predict with notable accuracy the stresses and strains that given loads induce in a patient's bones. We also show how these predictions can help to estimate the strength of the bones. Furthermore, we show how new numerical methods will enable a complete automation of this modelling process in the near future, thus opening the way to a much easier clinical deployment.

Then, from the middle level of interest (the organ), we shall report on the research results achieved on the upper dimensional (whole body) and lower dimensional scales (tissue, cell, molecule).

For the body model, we shall see how the neuromuscular indeterminacy problem remains substantially unsolved, and describe the new directions that are being explored.

For the tissue level, we report on recent improvements in the experimental methods that can characterize the bone tissue simultaneously in its morphology, biomechanics, biology and biochemistry, and how the combination of these different perspectives can open unexplored modelling opportunities. We shall also describe the Living Human Project (http://www.livinghuman.org/), an initiative aimed at collecting the huge amount of information required to produce a detailed physiome of the skeletal system.

Last, but not least, we shall briefly summarize the early results in bone mechanobiology, which might help in creating a complete predictive model of skeletal functional competence, and which also accounts for bone adaptation and its pharmacological, physiological and environmental determinants.

2. A subject-specific model of skeletal functional competence

Predicting subject-specific functional competence of a bone by developing and exploiting a bone model (geometric description at the organ level) involves three more distinct problems, originating at different dimensional scales, and converging onto the organ-level model (figure 1) as follows:

1. the prediction of the loads that daily life imposes on the bone (boundary conditions, at the body level);

   

   • Download figure
   • Open in new tab
   • Download PowerPoint

2. the prediction of how the bone tissue deforms under given mechanical stresses and the risk of fracture associated with such a stress/strain state (constitutive equation and failure criterion, at the tissue level); and

3. the prediction of the transformations of the mineralized extracellular matrix that are induced by the cellular activity, and their biomechanical determinants (bone remodelling, at the cell level).

3. An assessment protocol for predictive modelling technology

As quite a novel type of medical technology, subject-specific predictive modelling does not have an established assessment protocol yet. In industrial engineering, the problem of ensuring that predictive models are capable of fulfilling their intended purpose has been tackled in the context of quality management systems such as ISO 9000, which revolve around the concepts of verification and validation. For an extensive discussion on this topic, see Roache (1998).

The extension of these concepts to the peculiar context of computational biomechanics is not trivial, and found some attention in the specialized literature only recently (Viceconti et al. 2005; Anderson et al. 2007).

An editorial published in Clinical Biomechanics in 2005 (Viceconti et al. 2005) suggested extending the classic verification and validation approach to include two additional assessment stages before predictive models can be used in the clinical context: clinical accuracy and risk–benefit analysis. Clinical accuracy assessment is necessary to ensure that the accuracy of the predictive models observed in the ideal conditions of validation studies does not severely deteriorate in real-world clinical applications. Risk–benefit analysis was proposed as a means to evaluate if a new health care technology based on predictive models is, beyond its accuracy, truly useful in solving clinically relevant problems.

In the present paper, this paradigm is followed to report the state of the art in each sub-modelling activity, with only one major modification: consistently with the most recent health technology assessment approaches (e.g. NHS 2004), we refer to the last stage of the assessment protocol as clinical impact. Clinical impact includes risk–benefit analysis as a particular case. Thus, in the following sections, we shall evaluate the maturity of each sub-modelling activity according to these eight steps:

1. Observation. The description of the phenomenon to be described and/or predicted, and the quantitative measurements against which the model predictions shall be compared.

2. Model. A schematic description of the phenomenon under investigation, based on certain hypotheses that are needed to define the model. The model is often used to predict the phenomenon and therefore it is referred to as a ‘predictive model’.

3. Identification. The description of the input parameters of the model. If it is possible to measure an input parameter accurately, then it can be defined deterministically. Otherwise, if the parameter cannot be measured, statistical models are needed: a probabilistic description of the input parameters is established.

4. Verification. The procedure to make sure that the mathematical model is solved numerically with a known accuracy.

5. Sensitivity analysis. The identification and quantification of how the uncertainty affecting the input parameters propagates to the model predictions. The error affecting the measurement of some input parameters should not severely propagate to the model predictions. In addition, in all cases when other parameters cannot be quantified deterministically and need to be modelled as a true random factor (population indeterminacy), sensitivity should be determined.

6. Validation. Before a predictive model can be used to draw any conclusions, it is necessary to determine the accuracy with which the model is able to replicate the results of a controlled experiment. In order to attain a safe estimate of the model accuracy, the validation experiment should be designed so as to reproduce the conditions of intended use of the model (including the most challenging ones).

7. Clinical accuracy. The best validation experiment is designed in ideal conditions. However, in clinical practice, the true accuracy of the model (i.e. its capability to correctly diagnose positive and negative cases) can be much lower and needs to be assessed; in the context of diagnosis and prognosis, accuracy should be determined in terms of sensitivity and specificity (Everitt 1995).

8. Clinical impact. Even if it can be proved that subject-specific models are clinically accurate, this does not necessarily imply that they are useful. The ultimate criterion to assess the efficacy of a model is its clinical impact (e.g. if the use of a predictive model has the potential to change clinical practice, and if the prognosis for a given disease is significantly changed after the introduction of the predictive model).

The scope of this paper is not to provide a complete review of the existing literature on the issues above, as this would require far more space than a single manuscript. Instead, we aim to provide an overview of a state-of-the-art workflow for a specific application (mechanical integrity of the proximal human femur). To do so, the approaches taken by our group and by others are reported as an example.

4. Organ level: the bone model

The bone model is the core of the protocol for assessing and predicting functional competence, and is the starting point in this middle-out approach. The method of choice to model organ-level bone structures is finite-element (FE) analysis, as it allows an estimation of the state of stress/strain for virtually any geometry, material properties and loading conditions. Only very recently have alternative approaches been proposed, based on different numerical techniques, such as the kernel particle method (Liew et al. 2002) or the natural element method (Doblare et al. 2005). However, such alternative techniques are far from reaching the maturity of the FE methods.

In the case of the human femur, the core observation was that fractures are most frequent in the proximal metaphysis (Jeffery 1974; Rockwood et al. 1991; Rüedi & Murphy 2001). The position of the fracture can vary slightly, spanning the region from the femoral head (subcapital fracture) to the distal part of the trochanters (intertrochanteric fracture). They most frequently occur in elderly (osteoporotic) subjects, and are often associated with low-energy trauma (fall on the patient's side). A significant part (though not the majority) of the fractures of the untreated hip in the elderly are not associated with a primary traumatic event (Rockwood et al. 1991; Rüedi & Murphy 2001).

The main hypothesis concerning the proximal femur is that, in order to predict the risk of fracture for a given patient, one should estimate the force (or set of forces) that, when applied at a given location and in a given direction, does fracture the specific femur (Smith et al. 1992; Taylor 2003; Wang & Puram 2004; Schileo et al. 2008). Bone tissue is brittle in nature (see failure criterion, in §7). This implies that a bone fails in a certain region in a sudden way when the loads applied at the body level cause an excessive state of tensile stress/strain. This causes the fracture to initiate at a certain spot, and in a certain direction, which the organ-level model should be able to replicate.

A model of the femur (at the organ level) is first identified by its geometry, which can be obtained, also in vivo, from medical imaging datasets of the region of interest (van Rietbergen et al. 1995; Ota et al. 1999; Taylor 2003; Keyak et al. 2005; Taddei et al. 2007). The method of choice for the definition of the bone boundary is segmentation of computed tomography (CT) data, which provide a good contrast between mineralized bone tissue and surrounding soft tissues. When properly calibrated, CT data provide information also on the spatial distribution of the mineral density of the bone at a continuum level (Taddei et al. 2007). The spatial resolution of state-of-the-art CT in vivo is typically some tens of millimetres. The geometry of the femur extracted from the CT dataset was used as the input for the generation of the FE mesh to be used for prediction of the risk of fracture. The simplest implementation of the tissue constitutive model, to provide the spatial distribution at the organ level of the mechanical properties of the bone tissue, directly relates Young's modulus of the bone tissue to its mineral density. Automatic software was developed to map the non-homogeneous spatial distribution of the bone tissue material properties onto the mesh by performing a numerical integration of the calibrated CT density field over each FE volume (Taddei et al. 2004). This method is adopted by several independent studies (Perillo-Marcone et al. 2004; Gupta et al. 2006; Taylor 2006). Similar to the vast majority of other proposed approaches (e.g. Weinans et al. 2000; Keyak et al. 2001; Gupta et al. 2004; Barker et al. 2005; Bitsakos et al. 2005), this method provides locally isotropic material properties. Although it is well known that bone tissue is anisotropic, since currently available imaging technologies cannot provide information on the tissue-level anisotropy apart from some very peripheral regions, anisotropy was not considered.

Before being used to derive clinical indications, a model needs to be verified and validated (Keyak et al. 1993; Dalstra et al. 1995; Lengsfeld et al. 1998; Gomez-Benito et al. 2005; Viceconti et al. 2005). The term verification is commonly used to indicate the process ensuring that the numerical model accurately predicts the results of the theoretical model it is based on, which means assessing its numerical accuracy. Different aspects were analysed, resulting in the following guidelines. The first verification of the subject-specific modelling strategy, developed in the past few years, was based on post hoc error indicators (Zienkiewicz & Zhu 1992) related to mesh refinement. Such indicators suggest that, for modelling the human femur, an element size of 2 mm represents a good trade-off between numerical accuracy and computational weight (Viceconti et al. 2004; Decking et al. 2006; Gupta et al. 2006; Taylor 2006; Schileo et al. 2007). In addition, it was shown that the inclusion of the non-homogeneous material properties does not worsen the convergence rate of the model (Taddei et al. 2004).

Subsequently, the sensitivity of the model to the uncertainties in modelling the geometry and the material properties of a human femur from the CT data needs to be assessed (Keyak et al. 1993; Lengsfeld et al. 1998; Bayraktar et al. 2004; Taddei et al. 2006a). A sensitivity analysis, based on a Monte Carlo method, was performed on three femur models generated from the in vivo CT datasets (Taddei et al. 2006c). The geometry, the density and the mechanical properties of the bone tissue were considered as random input variables to take into account both measurement uncertainty and population variability. In fact, to assess the sensitivity of the FE models to such random input variables, the FE model outputs typically used in bone biomechanics (e.g. maximum principal strain, maximum principal stress, maximum strain energy density) can be treated as statistical variables (Pancanti et al. 2003; Belenguer Querol et al. 2006; Nicolella et al. 2006; Taddei et al. 2006c). The sensitivity of such statistical output variables to the variability of the inputs was assessed. The results (Taddei et al. 2006c) showed that, using the proposed method to build an FE model of a femur from a CT dataset of the quality typically achievable in clinical practice, the coefficients of variation in the output variables never exceed 9 per cent.

The following step was the validation of the modelling strategy against controlled experiments in the physiological range of loads (Keyak et al. 1993; Lengsfeld et al. 1998; Bayraktar et al. 2004; Yosibash et al. 2007a,b). Strain on the bone surface has successfully been measured by means of strain gauges (e.g. Lanyon et al. 1975; Huiskes et al. 1981; Rubin & Lanyon 1982; Field & Rushton 1989; Aamodt et al. 1997; Cristofolini 1997), as well as using optical techniques (e.g. Ferré et al. 1990; Katz et al. 1992; Spirakis et al. 1992; Gajda et al. 1995; Tyrer et al. 1995). To validate the FE modelling technique adopted by our group, a total of eight human femurs were tested in vitro (Cristofolini et al. 2006): bone strains were measured at a large number of points, while a set of loading scenarios that covered the range spanned in daily activities was applied (Bergmann 2001; Bergmann et al. 2001). Dedicated transducers were built to enable accurate control of the boundary conditions (Juszczyk et al. in press), to ensure that the applied loads in vitro replicated the external loading applied at the body level (see §5). The predicted principal strains were highly correlated with the experimentally measured ones (R²=0.91, regression line slope and intercept not significantly different from one and zero, respectively; Schileo et al. 2007). The r.m.s.e. was lower than 10 per cent of the maximum measured strain (Schileo et al. 2007). Similar results were reported by Yosibash et al. (2007a,b), but with a significantly lower number of experimental data for validation. When a failure criterion (see §7) based on the maximum principal strain (Taylor 2003; Bayraktar et al. 2004) was used to investigate femoral neck fracture replicated in vitro under physiologically relevant loading conditions, the model was able to identify correctly both the critical load and the location of the onset of failure (Schileo et al. 2008) under the conservative hypothesis that the whole bone would fracture once the yield strain was overcome at least in one location of the outer surface (figure 2). The accuracy of the bone superficial strain predictions was confirmed also in surgically treated femurs.

• Download figure
• Open in new tab
• Download PowerPoint

In order to assess the clinical accuracy of the organ-level model, one should theoretically perform a clinical study in which the patients enrolled in the study are examined, when still unfractured, using the proposed modelling technique (the probabilistic FE approach) to predict the actual risk of fracture under daily living activities (Cody et al. 1999; Davis et al. 2007). Such a study could be either retrospective or prospective. Typically, the risk threshold is increased progressively to form a receiver operating characteristic (ROC) curve (Metz et al. 1973), where the sensitivity and the specificity of the model's predictions are plotted. The greater the area under the ROC curve, the better the clinical accuracy of the model. In addition, the probability that the model predictions produce a false positive or false negative could be computed, and the clinical implications can be evaluated with risk analysis methods. However, prospective studies require a number of years, while the retrospective studies require a database of clinical observations including all the data necessary for subject-specific modelling (e.g. CT scan data of the region of interest), which is quite rare to find. Our group is working to establish such a clinical trial. In the meantime, we are using a probabilistic FE modelling approach to gain some understanding on how the model predictions can vary across a simulated population. The variation ranges and distributions of relevant patient conditions (bone quality, body weight, level of activity) have been characterized and included in the FE model formulation. The risk of fracture predicted by the probabilistic FE model over the simulated population was similar to that reported in epidemiological studies (Melton et al. 1986, 2005).

In some applications, the clinical impact of FE models has been assessed (Cody et al. 1999; Defranoux et al. 2005; Keyak et al. 2005; Melton et al. 2005). In our application, the clinical impact of the model of the femur was evaluated by implementing two distinct clinical applications: tumour and hip prostheses. In the first case, the evolution of the strength of a reconstructed femur during follow-up in a 10-year-old child was studied. The risk of fracture evolution was estimated starting from a personalized musculoskeletal model. The correlation of this risk of failure with the level of mineralization of the growing bone (bone remodelling, see §8) in the first 3 years after surgery was investigated (Taddei et al. 2006b). Another clinical application, in combination with experimental tests, regarded the design optimization of a proximal femoral epiphyseal replacement, with the aim of minimizing the risks associated with various biomechanical failure scenarios (Martelli et al. 2006), such as femoral neck fracture, proximal femoral bone resorption, progressive thinning of the cortical bone in the femoral neck. In this way, the suitability of the implant for a larger clinical indication could be evaluated. The fact that the model provided guidelines for the rehabilitation strategy (in the first case) and for the designer/manufacturer (in the second case) confirms the clinical relevance of the model proposed.

5. The middle-out approach: up and down

Bone fracture takes place at the organ level and at all the lower levels (including the tissue, cell and molecule levels). To focus on clinically relevant bone fractures (Rockwood et al. 1991; Rüedi & Murphy 2001), they are investigated at the organ level. After the organ-level model was fully developed and assessed, we started examining the higher and lower hierarchical levels (figure 1). The boundary conditions imposed on the organ-level model were defined at an upper scale, that of the whole-body level. The constitutive equation and the failure criterion were defined at the tissue level. Bone remodelling (i.e. the changes in bone tissue over time) was defined at the cell level.

As we shall see in the following sections, the development of these other models, and their coupling into a full multiscale model, is still a work in progress. We shall provide a description of the work that has been done, as well as of the challenges that remain open, for each dimensional scale.

6. Body level: bone loading conditions

The body-level model is the largest scale model deployed in this approach. This is used to determine the loading conditions that are caused at the organ level, in a given subject. This model must provide the lower level model with all the subject-specific boundary conditions (forces applied to the organ-level bone model).

The observations to be found at the body level are the loading conditions. This model must provide the lower level model with all the boundary conditions (forces transmitted between adjacent bone segments, forces applied by the muscles to the bone segment). At present, directly measuring the muscle forces exerted in vivo by the muscles during motion is an extremely challenging task. Indirect measurements have sometimes been performed in the past (Jacob et al. 1982), but the accuracy of such measurements is difficult to assess. The mechanical properties of the muscle–tendon complex have been measured in vivo by applying controlled loads (Onambele et al. 2006), but this does not provide a direct indication of the forces exerted during specific motor tasks. Measurement of the muscle forces using electromyography (EMG) provides only qualitative results (De Luca 1997). The muscle forces have sometimes been estimated using numerical simulations based on either motion analysis or ground reaction forces (e.g. Crowninshield et al. 1978; Patriarco et al. 1981; Davy & Audu 1987; Cheal et al. 1993; Fitzsimmons et al. 1996; Graham et al. 2008). A few experiments have been carried out directly measuring joint reaction forces, using diverse approaches and sometimes providing inconsistent results (Davy & Audu 1987; Davy et al. 1988; Hodge et al. 1989, Taylor et al. 1997; Bergmann 2001; Bergmann et al. 2001, 2004; Taylor & Walker 2001).

Load profiles recorded with a telemetric instrumented prosthesis were published by Bergmann (2001) and Bergmann et al. (2001). In addition, external boundary conditions are obtained from movement analysis and force platforms (Heller et al. 2001).

The underlying hypothesis for building a body-level model was that it could be used to predict the forces applied to a given organ, for specific motor tasks. With this aim, a musculoskeletal model of the lower limbs was developed that comprised 9 bone segments (pelvis, femurs, tibiae and fibulae, patellae and feet), 6 joints (hips, knee and ankles), 84 muscles and 108 skeletal landmarks (Montanari et al. 2006; figure 3). The muscles were represented with linear segments starting from the muscle origins to their insertion. Muscles with broad origins or insertion areas were represented with more than one segment. The musculoskeletal model was then registered in space and time with motion data collected from the same subject (Leardini et al. 2006). Running a static optimization algorithm on this model, we were able to predict muscle and joint forces acting on the patient's bone during a given motor task (Viceconti et al. 2006).

• Download figure
• Open in new tab
• Download PowerPoint

Predicted joint forces and muscle activation patterns were in good agreement with in vivo measurements with telemetric hip prostheses (Heller et al. 2001) and EMG recordings. However, a true subject-specific validation of the model will not be possible until biomechanical research solves the challenge of measuring the muscle forces directly in vivo and non-invasively. This would allow an iterative improvement loop to be established, where methods are refined until the subject-specific predictive accuracy is sufficiently good. However, even if this major challenge were solved, the modelling approach described above, which is the current state of the art, will still be capable of predicting only the optimal excitation patterns that produce the required movement/effort with minimal muscle and articular forces, given the optimization function adopted. Conversely, the problem is only partially solved for those paraphysiological conditions occurring when the neuromotor control is sub-optimal and the skeleton is overloaded. Thus, to provide the information needed for predicting the risk of neck fracture, the body-level model should ideally be able, given a set of inputs that characterize the patient's lifestyle, the neuromuscular control strategy and the musculoskeletal anatomy, to predict the spectrum of all possible loads (physiological or paraphysiological) the patient could exert on the bone during any motor task (figure 3).

The development of this type of probabilistic musculoskeletal model for a specific patient requires a radically different approach that includes (i) monitoring lifestyle with wearable technology, (ii) new imaging methods that make it possible to identify in a clinical setting all the subject-specific anatomo-functional quantities required for the model and (iii) quantitative neurological tests that make it possible to identify a probabilistic model of neuromuscular activation to the specific condition of each patient. Each of these steps poses methodological problems that are currently unsolved and are being addressed by various research groups throughout the world.

7. Tissue level: constitutive equation and failure criterion

The tissue-level model should provide the organ-level model with the equations describing the biomechanical behaviour of the bone tissue under load. This is obtained by defining a mathematical model (constitutive equation; figure 4) that relates the stresses induced by the external forces in the tissues with the deformation of the tissue itself, and then by identifying the constants of this model. Additionally, a criterion (failure criterion) must be identified to indicate the physical quantity that controls tissue fracture and the threshold leading to fracture. If the biomechanical behaviour of a tissue is described with sufficient accuracy by a linear isotropic elastic constitutive equation, two constants are sufficient to identify fully the stress–strain relationships in three dimensions (constitutive equations); additional constants are needed to identify the failure criterion. When the tissue behaviour requires the inclusion of anisotropy and/or nonlinearity to achieve the expected predictive accuracy, describing the constitutive equation and identifying its constants becomes a very challenging problem. In addition, tissues from different anatomical regions frequently exhibit significantly different properties, imposing a stochastic approach to the problem (Bevill et al. 2007). Similar problems are posed by the identification of the failure criterion. In fact, the correlation between bone structure (at the tissue level) and microdamage (in the form of linear microcracks and diffuse damage) has recently been investigated in relation to patient age (Diab & Vashishth 2007). Histological and histomorphometrical analyses were conducted to investigate the in vivo accumulation and localization of damage morphologies. Damage morphology showed no correlation with bone geometry parameters (organ level) and exhibited distinct preferences with bone microstructure (tissue level).

• Download figure
• Open in new tab
• Download PowerPoint

Currently, the organ-level model uses an isotropic heterogeneous linear elastic constitutive equation that is identified with the subject-specific tissue density distribution derived from properly calibrated CT data. The relationship between the tissue density and the elastic constants was derived empirically (Morgan et al. 2003) from experimental results obtained from the tissues of a number of subjects and from selected anatomical sites (figure 4), and then generalized to all subjects and to all sites. However, the validity of this generalization remains to be proved. In addition, the heterogeneous model provides a long-range anisotropy, but in some cases the short-range anisotropy of the tissue morphology is also relevant, and this is currently neglected.

There are two possible evolutions: the first one is to develop accurate tissue-level models of a large number of tissue biopsies taken from a representative population properly stratified for age, sex, pathology and anatomical sites. This will enable the gradual building of an atlas of the bone tissue constitutive equation defined statistically as a function of the population parameters.

The other possibility is related to the promise of modern imaging technologies to produce bone tissue images with a resolution (tens of micrometres) sufficient to determine the full three-dimensional tissue morphology. Currently, the radiation doses involved allow these technologies to be used only in the peripheral skeleton (e.g. the wrist). However, it can be expected that, within the next 5–10 years, imaging technology will be able to produce three-dimensional images of any skeletal region at tissue resolution with clinically acceptable radiation dose.

In both cases, tissue-level models need to be developed that are capable of predicting the tissue-level strain induced by a known boundary condition. The complexity of bone tissue morphology, especially the trabecular bone tissue (Gibson 2005), poses major challenges in creating FE models of portions of tissue. Current strategies involve the use of voxel-based mesh generators that cope with the topological complexity at the cost of a very large number (up to 100 million) of degrees of freedom to be solved (Ruimerman et al. 2005). This has pushed the development of special solvers capable of solving these problems in a reasonable amount of time using homogeneity assumptions (van Rietbergen et al. 1995; Niebur et al. 2000) or by exploiting massive parallelism (Arbenz et al. 2007). A more recent approach, developed by our group, uses a different numerical method to solve the elasticity equations, called the cells method (Taddei et al. in press). This method allows meshless implementations and exhibits a convergence rate much higher than that of the FE method. This enables implementations that couple the same level of automation as voxel-based mesh methods to a much smaller computational footprint for the model solution.

The second problem is the validation of these models. Some authors validate their model by comparing their structural displacement predictions with those measured on the same tissue biopsies used to generate the model. However, the cross-head displacement in a compression test is an extensive property of the structure, which can be obtained by infinite strain distributions. A true validation is possible only by measuring the displacement induced by a known load at a significant number of tissue points. High-resolution optical methods (Verhulp et al. 2004) can provide accurate measurement of the displacement field on the exposed surface of the specimen. Inspection of the internal structure of a tissue specimen with a spatial resolution of the order of 1 μm is possible by means of micro-computed tomography (micro-CT) imaging (Ruegsegger et al. 1996). The use of repeated micro-CT imaging during a compression test enables three-dimensional observation of the displacement field (Nazarian & Muller 2004). In practice, the development of a fully validated tissue-level model remains an open challenge.

A final aspect is the effect of the biochemical composition of the tissue. All tissue-level modelling studies assume that the properties of the tissue are constant and homogeneous. Actually, mineralized tissues are lamellar composites of collagen fibres organized in a complex three-dimensional fashion, stiffened by calcium phosphate crystals, which in part are bonded at the molecular level and in part are simply apposed onto the collagen fibre mesh. Significant alterations in the synthesis of the collagen, its spatial organization or the mineralization process may produce tissues that look morphologically equivalent but are functionally very different. But, as we said, this aspect has been neglected so far.

8. Cell level: bone remodelling

The tissue-level model translates the three-dimensional morphology of the bone tissue and its biochemical composition into a constitutive equation at the organ level. Even assuming that we are able to reconstruct the actual bone tissue morphology in a selected region of the patient's skeleton through high-resolution medical imaging, this morphology and the composition and arrangement of the tissue itself might change over time as a result of so-called bone remodelling.

All living tissues change their properties over time in relation to a number of factors (e.g. metabolism, pathologies, environment, exercise). This is caused by the interaction, at a sub-tissue level, of different cell types through a variety of molecules. Incorporating bone remodelling means incorporating a time-dependent model that updates the bone tissue morphology and its biochemical composition (and thus its constitutive equation) over time.

In the specific case of bone tissue, the resorption of existing extracellular matrix is performed by a multinuclear cell type called an osteoclast, formed by the fusion of cells of the monocyte–macrophage cell line, and characterized by high expression of tartrate-resistant acid phosphatase. The synthesis of collagen type I, which forms new bone tissue after the mineralization process, is performed by the osteoblasts, which arise from osteoprogenitor cells located in the periosteum and the bone marrow. Once osteoprogenitors start to differentiate into osteoblasts, they begin to express a range of proteins, including collagen type I, alkaline phosphatase, osteocalcin, osteopontin and osteonectin. In healthy conditions in the adult, the continuous activity of these two cellular types does not produce any net change in bone mass. However, metabolic changes related to ageing or the hormonal balance, alterations in environmental conditions including nutrition and lifestyle, or the emergence of pathological conditions may unbalance the process. Bone might respond to different stimuli such as variation in load and/or oxygenation by activation or depression of a cell metabolic pathway. These complex and often concurring mechanisms may result in an increase or a reduction in bone mass. In particular, in osteomalacia or osteoporosis, a progressive reduction in bone mass is observed.

This brief description is a crude simplification of the complex dynamical process that regulates bone metabolism. For example, osteoclasts and osteoblasts are known to interact with each other, for example, via the receptor activator of nuclear factor kappa B ligand. Similarly, the number of studies that attempt to model a selected portion of the bone mass regulation mechanism is too large to be detailed here. Consistently with the scope of this paper, attention here is on the approaches focusing on the mechanobiological aspects of the regulation mechanism.

Algorithms for interpreting and predicting bone remodelling at the tissue level have been published, based on both a simplistic black box approach and phenomenological observation (Cowin 1981; Cowin et al. 1994; Prendergast & Taylor 1994; Huiskes 1995; Luo et al. 1995; Taylor & Prendergast 1997; Huiskes et al. 1998; Fyhrie & Kimura 1999; Subbarayan & Bartel 2000).

A cell tissue-coupled model was proposed (Muller 2005), which predicted the changes in bone tissue morphology as a function of the osteoclastic penetration depth and apposition efficiency of the osteoblasts. The simulation itself describes an iterative process with a cellular remodelling cycle of 197 days (Muller 2005). If this model is coupled with the higher scale (the aforementioned tissue-level model), and with the lower scale (with models of cellular networks that predict how osteoclastic and osteoblastic activities are modulated by metabolic, pharmacological or environmental determinants), it could be possible to predict how alteration in such determinants affects the strength of the skeleton. In a recent interesting example of these cell-level models, it is suggested that Hill mechanisms of allosteric regulation could also be used in this context (Moroz & Wimpenny 2007).

9. Discussion

The aim of this paper was to discuss how a multiscale approach could be developed to investigate structural problems related to the human skeleton. This entire process is described by taking the strength (functional competence) of the proximal femur as an example. A middle-out approach was taken to describe the entire approach, so that the focus remains on the organ-level model, while integrating body-, tissue- and cell-level models.

In this brief overview, we showed how we established a fully validated subject-specific modelling strategy to develop organ-level models that can be used in clinical practice to predict the biomechanical competence of selected regions of the skeleton. Starting from this dimensional scale, we elaborated how the model could be extended up to the whole-body level and down to the tissue and cell levels. This constitutes a comprehensive predictive tool that can account also for other relevant factors such as alterations in neuromotor control, changes in tissue morphology or pharmacologically induced alterations in the bone remodelling process.

The resulting multiscale model (figure 5) is far from being completed. The farther from the organ-level scale, the more complex the problem becomes and the less advanced the state of the art is. However, we hope that this paper can provide a strategic vision that might guide the future work of those who are interested in this very exciting application of the physiome concept.

• Download figure
• Open in new tab
• Download PowerPoint