© Liò et al.; licensee BioMed Central Ltd. 2012
Published: 7 September 2012
Skip to main content
© Liò et al.; licensee BioMed Central Ltd. 2012
Published: 7 September 2012
This work focuses on the computational modelling of osteomyelitis, a bone pathology caused by bacteria infection (mostly Staphylococcus aureus). The infection alters the RANK/RANKL/OPG signalling dynamics that regulates osteoblasts and osteoclasts behaviour in bone remodelling, i.e. the resorption and mineralization activity. The infection rapidly leads to severe bone loss, necrosis of the affected portion, and it may even spread to other parts of the body. On the other hand, osteoporosis is not a bacterial infection but similarly is a defective bone pathology arising due to imbalances in the RANK/RANKL/OPG molecular pathway, and due to the progressive weakening of bone structure.
Since both osteoporosis and osteomyelitis cause loss of bone mass, we focused on comparing the dynamics of these diseases by means of computational models. Firstly, we performed meta-analysis on a gene expression data of normal, osteoporotic and osteomyelitis bone conditions. We mainly focused on RANKL/OPG signalling, the TNF and TNF receptor superfamilies and the NF-kB pathway. Using information from the gene expression data we estimated parameters for a novel model of osteoporosis and of osteomyelitis. Our models could be seen as a hybrid ODE and probabilistic verification modelling framework which aims at investigating the dynamics of the effects of the infection in bone remodelling. Finally we discuss different diagnostic estimators defined by formal verification techniques, in order to assess different bone pathologies (osteopenia, osteoporosis and osteomyelitis) in an effective way.
We present a modeling framework able to reproduce aspects of the different bone remodeling defective dynamics of osteomyelitis and osteoporosis. We report that the verification-based estimators are meaningful in the light of a feed forward between computational medicine and clinical bioinformatics.
There are two main types of bone tissues: cortical bone, and trabecular bone. The former is a compact tissue that makes up the outer shell of bones. It consists of a very hard (virtually solid) mass of bony tissue arranged in concentric layers called Haversian systems. Trabecular (also known as cancellous or "spongy") tissue is located beneath the compact bone and consists of a meshwork of bony bars (trabeculae) with many interconnecting spaces containing bone marrow. Both bone tissues undergo a continuous remodelling dynamics where old bone is replaced by new tissue ensuring the mechanical integrity and the morphology of the bone [1, 2]. However, pathological conditions such as cancer, infection and autoimmune diseases can alter the equilibrium between bone resorption and bone formation, reducing bone density and increasing the risk of spontaneous fractures.
Bone remodelling (BR) is a cellular process conducted by osteoclasts, the cells responsible for bone resorption and by osteoblasts, the cells responsible for bone formation. Osteoblasts follow osteoclasts in a highly coordinated manner, forming the so-called Basic Multi-cellular Units (BMUs). While osteoblasts and osteoclasts are located in the fluid part of the BMU, another type of cells, the osteocytes, are trapped in the bone matrix and they play a relevant role in the remodelling process. Osteocytes serve as mechanosensors: they translate mechanical stimuli at the tissue level into biochemical signals that flow through the osteocytic canalicular network to the BMU cells. In normal bone, the number of BMUs, the bone resorption rate, and the bone formation rate are all relatively constant .
Origination. During normal turnover or after a micro-crack, or as a response to mechanical stress, the osteocytes in the bone matrix produce biochemical signals showing sufferance towards the lining cells, i.e. the surface cells around the bone. The lining cells pull away from the bone matrix, forming a canopy which merges with the blood vessels.
Osteoclast recruitment. Stromal cells divide and differentiate into osteoblasts precursors. Pre-osteoblasts start to express RANKL, inducing the differentiation of and attracting pre-osteoclasts, which have RANK receptors on their surfaces. RANKL is a homotrimeric molecule displayed on the membrane of osteoblasts that stimulates differentiation in osteoclasts and is a key induction molecule involved in bone resorption leading to bone destruction.
Resorption. The pre-osteoclasts enlarge and fuse into mature osteoclasts. In cortical BMUs, osteoclasts excavate cylindrical tunnels in the predominant loading direction of the bone, while in trabecular bone they act at the bone surface, digging a trench rather than a tunnel. After the resorption process has terminated, osteoclasts undergo apoptosis.
Osteoblast recruitment. Pre-osteoblasts mature into osteoblasts and start producing osteoprotegerin (OPG). OPG inhibits the osteoclastic activity by binding to RANKL and preventing it from binding to RANK. When RANKL expression is high, osteoprotegerin levels are low and vice versa.
Mineralization. Osteoblasts fill the cavity by secreting layers of osteoids. Once the complete mineralization of the renewed tissue is reached, some osteoblasts can go apoptosis, other can turn into lining cells, while other can remain trapped in the bone matrix and become osteocytes.
Resting. Once the cavity has been filled by osteoblasts, the initial situation is re-established.
The bone remodelling undergoes a pathological process, generally related to ageing, termed osteopenia and with more severity, osteoporosis, during which an unbalance of the RANKL/OPG signalling equilibrium is typically observed. The osteoporosis is a skeletal disease characterized by low Bone Mineral Density (BMD) and structural fragility, which consequently leads to frequent micro-damages and spontaneous fractures; it is a chronic disease requiring long-term treatment. This disease primarily affects middle-aged women and elderly people and at present its social and economic impact is dramatically increasing, so much that the World Health Organization considers it to be the second-leading healthcare problem. While under normal circumstances, the ratio of RANKL/OPG is carefully balanced, the increase of RANKL plays an essential role in favouring resorption through osteoclast formation, function, and survival. With ageing and after a large number of remodelling cycles, the density of osteons increases and the cortical porosity and architectural defects of the bone increase as well. This leads to a vicious cycle where microdamages and consequently remodelling occur more and more frequently, weakening the bone structure and increasing the rate of spontaneous fractures . Moreover, recent studies suggest that plasma levels OPG and RANKL are inversely related to bone mineral density and contribute to the development of osteoporosis in postmenopausal women , and thalassemia-induced osteoporosis . One of the most worrying events is the infection of the bone which causes a disease called osteomyelitis. Similarly to osteoporosis, it is characterized by severe and rapid bone loss and by an unbalance at the molecular signalling level.
The aim of this work is to provide a computational modelling framework able to reproduce and compare the defective dynamics of osteoporosis and osteomyelitis. We believe that this framework could easily be adapted to model also other bone diseases like multiple myelomas or Paget's disease, and that it could help in better understanding the disruptions of cellular and signalling mechanisms that underlie such bone pathologies.
Osteomyelitis is a bone infection mainly caused by the aggressive pathogen S. aureus. Upon exposure to the bone, S. aureus induces a severe inflammatory response followed by progressive bone destruction and loss of the vasculature and with a persistent chronic infection; this is further complicated by the rapid emergence of resistant strains of S. aureus. Lab researches have shown that the infection prevents proliferation, induces apoptosis and inhibits mineralisation of cultured osteoblasts. The action of S. aureus increases RANKL expression and decreases OPG expression in osteoblasts in patients with staphylococcal osteomyelitis. Recent findings suggest that S. aureus SpA protein binds to osteoblasts, possibly through an interaction with the death receptor TNFR-1 which induces caspase 3 activation and apoptosis. The increase in RANKL is likely to trigger osteoclast-induced bone resorption and bone destruction and may help explain why patients with osteomyelitis have significant bone loss .
Although effective treatment of this disease is very difficult, one of most used drug is the fusidic acid that acts as a bacterial protein synthesis inhibitor by preventing the turnover of elongation factor G (EF-G) from the ribosome. Fusidic acid inhibits bacterial replication and does not kill the bacteria, and is therefore termed "bacteriostatic". Many strains of methicillin-resistant S. aureus (MRSA) remain sensitive to fusidic acid, but because there is a low genetic barrier to drug resistance (a single point mutation is all that is required), fusidic acid is usually combined with other antibiotics.
We believe that a model of the infection could provide a framework for a better diagnosis and understanding the antibiotic intervention. Here we develop a hybrid modelling framework for combining and untangling the relationships of physiological and molecular data. We then apply the methodology to determine disease related abnormalities of the key osteogenesis molecular network. The universality of the approach is demonstrated by an integration of the modelling and diagnosis which resembles medical visits with blood testing for infection progress and bone mineralisation measurements along a period of time. Our perspective is that this approach would inch towards an automatized methodology for improving disease classification and diagnosis.
Comparative representation of gene expression level for osteomyelitis and osteoporosis.
Regulation for Osteomyelitis (GPL96)
Regulation for Osteomyelitis (GPL97)
Regulation for Osteoporosis (GPL96)
Interestingly we found that, despite a very small increase of RANKL gene expression in osteoporosis and a larger increase in osteomyelitis, OPG gene expression become more deregulated in both osteomyelitis and osteoporosis. There is the increased expression of different isoforms of OPG which are known to have different binding capability with RANKL and seem to be linked, from mice experiments, to hypocalcemia . Therefore we report that gene expression in osteoporosis and osteomyelitis could generate an unbalance between RANKL and OPG due to the different OPG isoforms, but also other genes, related to TNF, TNF receptor superfamilies and to NF-kB may be involved. Although gene expression and actual protein abundance are only loosely correlated, taking into account the results of gene expression data, we modified the autocrine and paracrine parameters of the existing mathematical model based on Komarova model . We considered more appropriate to incorporate into the model the algebraic relationship of positive and negative regulators (such as RANKL and OPG) than just the RANKL change. On the basis of this consideration we developed new models for reproducing osteoporotic and osteomyelitis conditions.
In this work we present a combined computational framework for the modelling, simulation and verification of the bone remodelling process, and of bone pathologies like osteomyelitis and osteoporosis. Based on the methods developed in [11, 12], this approach consists of the following two building blocks:
We develop a differential equation model for describing the dynamics of bone remodelling and of bone-related pathologies at a multicellular level. The model describes the continuous changes of, and the interactions between populations of osteoclasts and osteoblast (including bacteria in the osteomyelitis model). Bone density is calculated as the difference between the formation activity which is proportional to osteoblasts concentration, and the resorption activity which is proportional to osteoclasts concentration. In the last twenty years, a variety of mathematical and computational models has been proposed in order to better understand the dynamics of bone remodelling (reviewed in [13–15]). Three main categories of models can be distinguished: those focusing on the organ level, where bone is described as a continuum material only characterized by its density; on the biomechanical properties and on the microstructural information at the tissue level; and on the cellular level where the interactions occurring among the different types of bone cells are concerned. The latter category can also incorporate intracellular signalling pathways and mechanosensing mechanisms (i.e. the process by which mechanical stimuli are translated into cellular signals). Our cellular-level model is based on the work by Komarova et al , where they developed an important model for BR based on experimental results described in Parfitt's work  which has inspired many other similar models. In particular we extended it in order to explicitly simulate bone pathologies: osteoporosis is reproduced by including an ageing factor that decreases the death rates of cells and by including a factor that increases the RANKL expression; osteomyelitis is modelled by adding a state variable for bacteria that affects the autocrine and paracrine regulation factors of osteoblasts and osteoclasts, similarly to Ayati's model on bone myeloma . Although several efforts have been made in developing mathematical model for osteomyelitis and osteoporosis, molecular data has been rarely considered so far, regardless the availability of different gene expression microarray data related to osteomyelitis and osteoporosis and based on only single microarray database. So, we have developed mathematical model and showed the comparative study of gene expression data from different databases of similar platform to find out the genes expression level related to the RANKL, RANK, OPG and NF-kB proteins, which are strongly related to the osteomyelitis and osteoporosis.
We define a stochastic model for bone remodelling from the ODE specification, that allows us to analyse the random fluctuations and the discrete changes of bone density and bone cells. Given that randomness is an inherent feature of biological systems, whose components are naturally discrete, the stochastic approach could give useful insights on the bone remodelling process. Indeed, stochasticity plays a key role in bone remodelling, e.g. the fluctuations in molecular concentrations of RANKL and OPG produce changes in the chemotaxis (the process by which cells move toward attractant molecules) of osteoclasts and osteoblasts. This may affect for example the cell differentiation, number and arrival time, and consequently the whole remodelling process. Besides achieving a good fitting between the ODE model and the stochastic one, we employ probabilistic model checking techniques for deriving three different clinical estimators that enable to assess the expected bone density, the density change rate, and the variance of bone density. Model checking is a static technique for automatically search for a property (specified as a logical formula) to hold or not over a definite set of states, and relies on qualitative properties: given a model and a property to verify, it returns an affirmative or a negative answer, i.e. the property holds or not. Differently, probabilistic model checking is equipped with quantitative information, and given a stochastic model and a property to verify, it returns the probability of the formula being satisfied. We believe that this kind of quantitative, formal and automated analysis may represent a step ahead in the understanding of bone diseases like osteomyelitis and osteoporosis, by shifting the attention from an informative, but empirical, analysis of the graphs produced by simulations towards more precise quantitative interpretations.
The model describes the autocrine and paracrine relationships between osteoclasts and osteoblasts. Autocrine signalling usually occurs by a secreted chemical interacting with receptors on the surface of the same cell. In the paracrine process a chemical signals that diffuse outside the emitting cell and interacts with receptors on nearby cells. Here the parameters g ij describe the effectiveness of autocrine and paracrine regulation, s.t. g 11 describes the osteoclast autocrine regulation, g 22 the osteoblast autocrine regulation, g 21 is the osteoblast-derived paracrine regulation, and g 12 is the osteoclast paracrine regulation. The nonlinearities of these equations are approximations for the interactions of the osteoclast and osteoblast populations in the proliferation terms of the equations. The autocrine signalling has a positive feedback on osteoclast production (g 11 >0), and paracrine signalling has a negative feedback on osteoclast production (g 21 <0). The autocrine signalling has a positive feedback on osteoblast production (g 22 >0), and paracrine signalling has a positive feedback on osteoblast production (g 12 >0).
where and denote the steady states of O c and O b , resp. For the spongy type bone we consider the variable z as the localized trabecular mass beneath a point on the bone surface.
In order to reproduce the defective dynamics (i.e. bone negative balance) characterizing osteoporosis, we assumed an increased death rate for osteoclasts and osteoblasts, motivated by the fact that the occurrence of defective bone pathologies in elderly patients is partly attributable to the reduced cellular activity typical of those patients. Therefore we introduced the parameter gageing as a factor multiplying the death rates β i .
This model has been inspired from Ayati's work on multiple myeloma bone disease  and the key difference with respect to Komarova's model  is the addition of the terms f ij B/s that couple the bacterial density and its maximum size to the power laws for the osteoclast/osteoblast interactions. The bacterial parameters f 11, f 12, f 21, f 22 are all nonnegative. The S. aureus-induced infection affects the normal remodelling activity by:
reducing osteoblasts' growth rate: in fact, the paracrine promotion of osteoblasts is reduced , and the autocrine promotion of osteoblasts is reduced as well ;
increasing RANKL and decreasing OPG expression: as previously stated, the paracrine inhibition of osteoclasts is a negative exponent resulting from the difference between the effectiveness of OPG signalling and that of RANKL signalling. Since , the infection causes an increase in RANKL expression and therefore a decrease in OPG expression.
In addition the infection increases the autocrine promotion of osteoclasts (since g 11 >0). We have taken γ B to be independent of bone loss. The parameter V describes the effectiveness of the antibiotic treatment as a factor decreasing the growth rate γ B of bacteria. Two different kinds of treatment can be distinguished: bacteriostatic treatments that stop bacteria proliferation (V = γ B ); and bacteriocide treatments which kill bacteria (V > γ B ).
O c and O b growth rates
(3, 4) day -1
O c and O b death rates
(0.2, 0.02) day -1
(g 11, g 12, g 22, g 21)
Effectiveness of autocrine/paracrine regulation
(1.1, 1, 0, -0.5)
(k 1, k 2)
Resorption and formation rates
(0.0748, 0.0006395) day -1
S. aureus growth rate
0.005 day -1
S. aureus carrying capacity
Effectiveness of antibiotic treatment
(0.005, 0.007) day -1
(200, 400, 600) days
(f 11, f 12, f 22, f 21)
Effect of infection on regulation factors
(0.005, 0, 0.2, 0.005) day -1
Steady levels of O c and O b
Control: (1.16, 231.72)
Osteoporosis: (1.78, 177.91)
Osteomyelitis: (5, 316)
(O c0, O b0, B 0)
Control: (11.16, 231.72, 1)
Osteoporosis: (11.78, 177.91, 1)
Osteomyelitis: (15, 316, 1)
Following and extending the work in , we define a stochastic model for bone remodelling and perform formal analysis by means of probabilistic verification techniques, which allow to assess the probability of a particular configuration of the biological system (usually expressed as a logical formula) being reached. In our settings we derive a Continuous Time Markov Chain (CTMC) from the mathematical model described above and we use the model checker PRISM .
where x max and x min are the maximum and the minimum x, resp. In other words, growth rates in the ODE model become the stochastic rates in a transition incrementing the population size, while death rates are involved in the transitions decrementing the population size.
Transitions in the stochastic model for bone remodelling.
: O c = O c + 1
O c > 0 →
gageingβ1 O c
: O c = O c - 1
O c > 0 →
k 1 O c
: O b = O b + 1
O b > 0 →
gageing β2 O b
: O b = O b - 1
O b > 0 →
k 2 O b
0 <B <B max ∧ treat = 0 →
: B = B + 1
treat = 0 →
: treat = 1
0 <B <B max ∧ treat = 1 ∧ V < γ B →
: B = B + 1
B > 0 ∧ treat = 1 ∧ V > γ B →
: B = B - 1
(d) Bone resorbed reward
(e) Bone formed reward
[resorb] true: 1
[form] true: 1
Osteomyelitis and osteoporosis are assessed through the verification of quantitative properties over the defined stochastic model.
Let assume that the simulation of the PRISM implementation of the model is run in parallel with the determination of clinical parameters during the periodic medical visits of a patient. These medical visits provide a mean of fine tuning a personalised model of the disease and a measure of how a therapy is effective. Different diseases, when monitored in a continuous way, may produce different alterations in local mineral density. We could extend the statistical estimators of a disease to: 1) the BMD (measured as z-score, the number of standard deviations above or below the mean for the patients age, sex and ethnicity; or as t-score, i.e. the number of standard deviations above or below the mean for a healthy 30 year old adult of the same sex and ethnicity as the patient); 2) The rate of change of BMD. This estimator tells us the emergence of defects of the bone metabolism in terms of signaling networks of RANK/RANKL and decrease of pre-osteoblast number; 3) The variance, skewness and curtosis of the the local small scale intermittency of the signal. For example osteomyelitis and osteoporosis show slightly confounding pattern of BMD decrease; we could also think at the confounding patterns of IRIS in HAART therapy, co-morbidity of osteopetrosis and osteoporosis, multiple myelomas, breast cancer, diabetes and metabolic syndromes, etc. The variance could perhaps help in discriminating among bone-related diseases. From a technical viewpoint, properties to verify have been formulated in CSL (Continuous Stochastic Logic) , and they give rise to three clinical estimators that we evaluate over 1200 days (about four years), which is enough to assess the presence of bone diseases:
Density variance. While the first estimator computes the expected value of bone density, here we calculate the variance of BMD taking into account the whole state space and the actual bone density at each state.
Here we report that the genetic complexity and the gene expression data meta analysis shows that the underlying "mystery" of bone remodelling is much greater than handled by the current mathematical models. In other words we are not able to use all our gene expression results in a full model of BR diseases. Although our model of osteomyelitis and the comparison with the osteoporosis is not able to consider all this complexity, nevertheless it makes a partial use of the results of the analysis of the experimental data and produces a realistic description of the pathology. From a methodological point of view the combination of mathematical and formal method approach has led to the proposal of considering additional estimators (first derivatives and variance) of the bone pathologies as diagnostic tool. That could also inspire the ideal situation in which a personalised model is generated from (personalised) data and the comparison between clinical data obtained during periodic medical check-up is compared with the computer predictions.
We found that there are no comprehensive analysis on osteomyelitis; most studies focus on specific conditions. We have collected a large ensemble of gene expression data related to osteomyelitis and osteoporosis. For this reason, we have considered 6 microarray data sets of the same platform GPL96 from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/), accession numbers are GSE16129, GSE6269, GSE11907, GSE11908, GSE13850 and GSE7429 [20–23]. We observe that RANKL, RANK, OPG and NF-kB proteins impact more on the bone remodelling for osteomyelitis and osteoporosis [7, 20–22]. For this reason to understand the effect osteomyelitis and osteoporosis on bone remodelling, we have considered the genes related to the proteins RANKL, RANK, OPG, NF-kB proteins, TNF and TNF receptor superfamilies. We observed that there are 82 genes are related with these proteins. So, we filtered the required 82 genes related data. We have selected samples for 48 infected and 27 healthy controls for osteomyelitis and 30 infected and 30 healthy controls for osteoporosis. The datasets contain data from people of different age and sex.
For more evidence about osteomyelitis, we have considered more gene expression data related to osteomyelitis on different platform GPL97. For this reason, we have considered additional 3 microarray data sets from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/), accession numbers are GSE6269, GSE11907 and GSE11908 [21, 22]. To understand the effect of osteomyelitis on the bone remodelling, we have considered the genes related to the proteins RANKL, RANK, OPG, NF-kB proteins, TNF and TNF receptor superfamilies like previous analysis. We observed that in the platform GPL97, there are 31 genes are related to these proteins and superfamilies. So, we filtered the required genes related data. We have selected samples for 43 infected and 17 healthy controls. Standard anova and Box plots representation were used to analyse and visualise the expression levels of these genes for the infection of osteomyelitis and osteoporosis condition. We output in Table 1 the groups of over expressed and under expressed categories.
We have implemented the ODE model based on Komarova et al  in R, and using the FME package  to analyse parameter sensitivity and robustness. We have used Mathematica and MATLAB for steady states and ODE calculation using state of art numerical routines. Scripts and functions for the models could be made available upon request to the first author. For the specification of the stochastic model and for performing probabilistic verification we have adopted the open-source PRISM probabilistic model checker , one of the reference existing model checkers for the analysis of systems which exhibits random or probabilistic behaviour. Since model checking is based on graph-theoretical techniques for exploring the whole state space of the model, this task becomes computationally infeasible for non-trivial models, due to the combinatorial explosion of the state space. For this reason, verification has been performed by means of approximate probabilistic model checking techniques that calculate the probability of a given property on a statistical basis, i.e. by sampling on a number of simulations of the model. In this work we have taken 20 samples for each verified property, that were enough to reproduce outputs similar to the non-approximate verification. PRISM models could be made available upon request to the second author.
We thank Bruce P. Ayati (Iowa University) and Glenn Webb (Vanderbilt University) for suggestions and help in computation.
This article has been published as part of BMC Bioinformatics Volume 13 Supplement 14, 2012: Selected articles from Research from the Eleventh International Workshop on Network Tools and Applications in Biology (NETTAB 2011). The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2105/13/S14
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.