Modelling osteomyelitis
 Pietro Liò^{1},
 Nicola Paoletti^{2}Email author,
 Mohammad Ali Moni^{1},
 Kathryn Atwell^{1},
 Emanuela Merelli^{2} and
 Marco Viceconti^{3}
https://doi.org/10.1186/1471210513S14S12
© Liò et al.; licensee BioMed Central Ltd. 2012
Published: 7 September 2012
Abstract
Background
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.
Results
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 metaanalysis 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 NFk B 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.
Conclusions
We present a modeling framework able to reproduce aspects of the different bone remodeling defective dynamics of osteomyelitis and osteoporosis. We report that the verificationbased estimators are meaningful in the light of a feed forward between computational medicine and clinical bioinformatics.
Keywords
Background
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 socalled Basic Multicellular 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 [3].
 1.
Origination. During normal turnover or after a microcrack, 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.
 2.
Osteoclast recruitment. Stromal cells divide and differentiate into osteoblasts precursors. Preosteoblasts start to express RANKL, inducing the differentiation of and attracting preosteoclasts, 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.
 3.
Resorption. The preosteoclasts 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.
 4.
Osteoblast recruitment. Preosteoblasts 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.
 5.
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.
 6.
Resting. Once the cavity has been filled by osteoblasts, the initial situation is reestablished.
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 microdamages and spontaneous fractures; it is a chronic disease requiring longterm treatment. This disease primarily affects middleaged 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 secondleading 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 [4]. 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 [5], and thalassemiainduced osteoporosis [6]. 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
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 TNFR1 which induces caspase 3 activation and apoptosis. The increase in RANKL is likely to trigger osteoclastinduced bone resorption and bone destruction and may help explain why patients with osteomyelitis have significant bone loss [7].
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 (EFG) from the ribosome. Fusidic acid inhibits bacterial replication and does not kill the bacteria, and is therefore termed "bacteriostatic". Many strains of methicillinresistant 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.
Results and discussion
Meta analysis of gene expression data
Comparative representation of gene expression level for osteomyelitis and osteoporosis.
Regulation for Osteomyelitis (GPL96)  Gene ID  Regulation for Osteomyelitis (GPL97)  Gene ID  Regulation for Osteoporosis (GPL96)  Gene ID 

Up regulated  NFKB2_1  Up regulated  NFKB2_1  Up regulated  NFKB1 
NFKB2_2  NFKBIZ_1  NFKB2_1  
NFKBIA  NFKBIZ_2  NFKB2_2  
NFKBIE  RELL1  REL_2  
REL_2  RELT  RELA_1  
RELB  TNFSF13B_1  RELA_2  
TNFRSF10B_2  TNFSF13B_2  RELB  
TNFRSF10C_2  TRAF7_1  TNFRSF17  
TNFRSF10C _3  TRAF7_3  TNFSF10_2  
TNFRSF10C_4  TRAFD1_2  TRAF3_1  
TNFRSF1A  Down regulated  TNFRSF10A  Down regulated  TNFRSF10B_2  
TNFRSF1B  TNFRSF18_2  TNFRSF25_2  
TNFSF10_1  TRAF1  TNFSF10_3  
TNFSF10_2  TRAF3IP1  TRAF3IP3_1  
TNFSF10_3  TRAF3IP3_1  TRAF3IP3_3  
TNFSF12_3  TRAF3IP3_2  TRAF5  
TNFSF12_4  
TNFSF12_2  
TNFSF13  
TRAF3IP3_2  
TRAF3IP3_3  
TRAFD1_2  
Down regulated  IKBKG 2  
NFKB1  
RELA_1  
TNFRSF14  
TNFRSF25_1  
TNFRSF25_2  
TNFRSF25_3  
TNFRSF25_4  
TNFRSF25_6  
TRAF1  
TRAF3IP2_2  
TRAF3IP3_1  
TRAF5 
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 [9]. 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 NFkB 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 [10]. 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.
A computational framework for bone dynamics
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:
Mathematical model
We develop a differential equation model for describing the dynamics of bone remodelling and of bonerelated 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 cellularlevel model is based on the work by Komarova et al [10], where they developed an important model for BR based on experimental results described in Parfitt's work [8] 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 [16]. 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 NFkB proteins, which are strongly related to the osteomyelitis and osteoporosis.
Model verification
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.
Modelling bone remodelling pathologies
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 osteoblastderived 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 ${\mathit{\u014c}}_{c}$ and ${\mathit{\u014c}}_{b}$ 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 g_{ageing} as a factor multiplying the death rates β_{ i }.
Osteomyelitis effects on bone remodelling
This model has been inspired from Ayati's work on multiple myeloma bone disease [16] and the key difference with respect to Komarova's model [10] 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. aureusinduced infection affects the normal remodelling activity by:

reducing osteoblasts' growth rate: in fact, the paracrine promotion of osteoblasts is reduced $\left({g}_{12}/\left(1+{f}_{12}\frac{B}{s}\right)<{g}_{12},\phantom{\rule{0.3em}{0ex}}\mathsf{\text{since}}\phantom{\rule{0.3em}{0ex}}{g}_{12}>0\right)$, and the autocrine promotion of osteoblasts is reduced as well $\left({g}_{22}{f}_{22}\frac{B}{s}<{g}_{22}\right)$;

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 ${g}_{21}\left(1{f}_{21}\frac{B}{s}\right)>{g}_{21}$, 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 }).
Model parameters.
Parameter  Description  Value 

(α_{1}, α_{2})  O_{ c } and O_{ b } growth rates  (3, 4) day^{1} 
(β_{1}, β_{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} 
g_{ageing}  Ageing factor  2 
g_{por}  RANKL factor  0.1 
γ_{ B }  S. aureus growth rate  0.005 day^{1} 
s  S. aureus carrying capacity  100 
V  Effectiveness of antibiotic treatment  (0.005, 0.007) day^{1} 
t _{ treat }  Dosage time  (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} 
$\left(\overline{{O}_{c}},\overline{{O}_{b}}\right)$  Steady levels of O_{ c } and O_{ b }  Control: (1.16, 231.72) 
Osteoporosis: (1.78, 177.91)  
Osteomyelitis: (5, 316)  
(O_{c 0}, O_{b 0}, B_{0})  Initial states  Control: (11.16, 231.72, 1) 
Osteoporosis: (11.78, 177.91, 1)  
Osteomyelitis: (15, 316, 1) 
Stochastic model for the verification of bone pathologies
Following and extending the work in [11], 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 [17].
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.
(a) Osteoclasts  

[]  $0<{O}_{c}<{O}_{{c}_{max}}\Lambda \phantom{\rule{0.3em}{0ex}}{O}_{b}>0\to $  ${\alpha}_{1}{O}_{c}^{{g}_{11}\left(1+{f}_{11}\frac{B}{s}\right)}{O}_{b}^{{g}_{21}\left(1{f}_{21}\frac{B}{s}\right)}$  : O_{ c } = O_{ c } + 1 
[]  O_{ c } > 0 →  g_{ageing}β_{1}O_{ c }  : O_{ c } = O_{ c }  1 
[resorb]  O_{ c } > 0 →  k _{1} O _{ c }  : true 
(b) Osteoblasts  
[]  $0<{O}_{b}<{O}_{{b}_{max}}\Lambda \phantom{\rule{0.3em}{0ex}}{O}_{c}>0\to $  ${\alpha}_{2}{O}_{c}^{{g}_{12}\left(1+{f}_{11}\frac{B}{s}\right)}{O}_{b}^{{g}_{22}{f}_{22}\frac{B}{s})}$  : O_{ b } = O_{ b } + 1 
[]  O_{ b } > 0 →  g_{ageing} β_{2}O_{ b }  : O_{ b } = O_{ b }  1 
[form]  O_{ b } > 0 →  k _{2} O _{ b }  : true 
(c) Bacteria  
[]  0 <B <B_{ max } ∧ treat = 0 →  ${\gamma}_{B}B\cdot log\left(\frac{s}{B}\right)$  : B = B + 1 
[]  treat = 0 →  $\frac{1}{{t}_{treat}}$  : treat = 1 
[]  0 <B <B_{ max } ∧ treat = 1 ∧ V < γ_{ B } →  $\left({\gamma}_{B}V\right)B\cdot log\left(\frac{s}{B}\right)$  : B = B + 1 
[]  B > 0 ∧ treat = 1 ∧ V > γ_{ B } →  $\left(V\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{\gamma}_{B}\right)B\cdot log\left(\frac{s}{B}\right)$  : B = B  1 
(d) Bone resorbed reward  (e) Bone formed reward  
[resorb] true: 1  [form] true: 1 
Potentialities in clinical bioinformatics and conclusions
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 zscore, the number of standard deviations above or below the mean for the patients age, sex and ethnicity; or as tscore, 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 preosteoblast 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, comorbidity of osteopetrosis and osteoporosis, multiple myelomas, breast cancer, diabetes and metabolic syndromes, etc. The variance could perhaps help in discriminating among bonerelated diseases. From a technical viewpoint, properties to verify have been formulated in CSL (Continuous Stochastic Logic) [19], 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:

Bone density estimator. It is calculated as the difference between the cumulative $\left({\mathcal{C}}^{\le t}\right)$ expected rewards (ℛ_{{"..."}}) for bone formation and bone resorption, with the formula${f}_{BD}\left(t\right):{\mathcal{R}}_{\left\{\prime \prime \mathsf{\text{bone}}\phantom{\rule{0.3em}{0ex}}\mathsf{\text{Formed}}\prime \prime \right\}}=?\left[{\mathcal{C}}^{\le t}\right]{\mathcal{R}}_{\left\{\prime \prime \mathsf{\text{bone}}\phantom{\rule{0.3em}{0ex}}\mathsf{\text{Resorbed}}\prime \prime \right\}}=?\left[{\mathcal{C}}^{\le t}\right],\phantom{\rule{0.3em}{0ex}}t=0,\phantom{\rule{0.3em}{0ex}}10,...,1200$

Density change rate. It allows to assess rapid negative and positive changes in bone density. This estimator could be particularly helpful in detecting the insurgence of osteomyelitis before critical values of bone density are reached, since osteomyelitis is typically characterized by a higher negative change rate than osteoporosis. In particular the estimator is defined as the difference quotient of BMD over a time interval of months, e.g. 50 days. The formula obtained is$\frac{{f}_{BD}\left(t+\Delta t\right){f}_{BD}\left(t\right)}{\Delta t},\phantom{\rule{0.3em}{0ex}}t=0,10,...,1200.$

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 checkup is compared with the computer predictions.
Methods
Data analysis
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 NFkB 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, NFkB 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, NFkB 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.
ODE and probabilistic model checking models
We have implemented the ODE model based on Komarova et al [10] in R, and using the FME package [24] 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 opensource PRISM probabilistic model checker [17], one of the reference existing model checkers for the analysis of systems which exhibits random or probabilistic behaviour. Since model checking is based on graphtheoretical techniques for exploring the whole state space of the model, this task becomes computationally infeasible for nontrivial 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 nonapproximate verification. PRISM models could be made available upon request to the second author.
Declarations
Acknowledgements
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/bmcbioinformatics/supplements/13/S14
Authors’ Affiliations
References
 Manolagas S, Parfitt A: What old means to bone. Trends in Endocrinology & Metabolism. 2010, 21 (6): 369374. 10.1016/j.tem.2010.01.010.View ArticleGoogle Scholar
 Karsenty G, Oury F: The central regulation of bone mass, the first link between bone remodeling and energy metabolism. Journal of Clinical Endocrinology & Metabolism. 2010, 95 (11): 479510.1210/jc.20101030.View ArticleGoogle Scholar
 Raggatt L, Partridge N: Cellular and molecular mechanisms of bone remodeling. Journal of Biological Chemistry. 2010, 285 (33): 2510310.1074/jbc.R109.041087.PubMed CentralView ArticlePubMedGoogle Scholar
 Whitfield J: Growing bone. Landes Bioscience. 2007Google Scholar
 Jabbar S, Drury J, Fordham J, Datta H, Francis R, Tuck S: Osteoprotegerin, RANKL and bone turnover in postmenopausal osteoporosis. Journal of Clinical Pathology. 2011, 64 (4): 35410.1136/jcp.2010.086595.View ArticlePubMedGoogle Scholar
 Nea Morabito: Osteoprotegerin and RANKL in the Pathogenesis of ThalassemiaInduced Osteoporosis: New Pieces of the Puzzle. Journal of Bone and Mineral Research. 2004, 19 (5): 722727. 10.1359/jbmr.040113.View ArticleGoogle Scholar
 Claro T, Widaa A, O'Seaghdha M, Miajlovic H, Foster T, O'Brien F, Kerrigan S: Staphylococcus aureus Protein A Binds to Osteoblasts and Triggers Signals That Weaken Bone in Osteomyelitis. PloS one. 2011, 6 (4): e1874810.1371/journal.pone.0018748.PubMed CentralView ArticlePubMedGoogle Scholar
 Parfitt A: Osteonal and hemiosteonal remodeling: The spatial and temporal framework for signal traffic in adult human bone. Journal of cellular biochemistry. 1994, 55 (3): 273286. 10.1002/jcb.240550303.View ArticlePubMedGoogle Scholar
 He Z, Yang G, Chen Z, Li B, Zhang W, Wu X: A novel isoform of osteoprotegerin gene: cloning and expression and its hypocalcemic effect in mice. Protein and Peptide Letters. 2000, 7 (4): 233240.Google Scholar
 Komarova S, Smith R, Dixon S, Sims S, Wahl L: Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone remodeling. Bone. 2003, 33 (2): 206215. 10.1016/S87563282(03)001571.View ArticlePubMedGoogle Scholar
 Liò P, Merelli E, Paoletti N: Multiple verification in computational modeling of bone pathologies. CompMod. 2011, 8296.Google Scholar
 Paoletti N, Lio P, Merelli E, Viceconti M: Multilevel Computational Modeling and Quantitative Analysis of Bone Remodeling. IEEE/ACM Transactions on Computational Biology and Bioinformatics. 2012, 99 (PrePrints):Google Scholar
 Geris L, Vander Sloten J, Van Oosterwyck H: In silico biology of bone modelling and remodelling: regeneration. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2009, 367 (1895): 203110.1098/rsta.2008.0293.View ArticleGoogle Scholar
 Gerhard F, Webster D, van Lenthe G, Müller R: In silico biology of bone modelling and remodelling: adaptation. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2009, 367 (1895): 201110.1098/rsta.2008.0297.View ArticleGoogle Scholar
 Pivonka P, Komarova S: Mathematical modeling in bone biology: From intracellular signaling to tissue mechanics. Bone. 2010, 47 (2): 181189. 10.1016/j.bone.2010.04.601.View ArticlePubMedGoogle Scholar
 Ayati B, Edwards C, Webb G, Wikswo J: A mathematical model of bone remodeling dynamics for normal bone cell populations and myeloma bone disease. Biology Direct. 2010, 5: 2810.1186/17456150528. [http://www.biologydirect.com/content/5/1/28]PubMed CentralView ArticlePubMedGoogle Scholar
 Kwiatkowska M, Norman G, Parker D: PRISM 4.0: Verification of Probabilistic Realtime Systems. Proc 23rd International Conference on Computer Aided Verification (CAV'11), Volume 6806 of LNCS, Springer. 2011, 585591.Google Scholar
 Dayar T, Mikeev L, Wolf V: On the numerical analysis of stochastic LotkaVolterra models. Computer Science and Information Technology (IMCSIT), Proceedings of the 2010 International Multiconference on, IEEE. 2010, 289296.Google Scholar
 Aziz A, Sanwal K, Singhal V, Brayton R: Model checking continuous time Markov chains. ACM Trans Computational Logic. 2000, 1: 162170. 10.1145/343369.343402.View ArticleGoogle Scholar
 Ardura M, Banchereau R, Mejias A, Di Pucchio T, Glaser C, Allantaz F, Pascual V, Banchereau J, Chaussabel D, Ramilo O: Enhanced monocyte response and decreased central memory T cells in children with invasive Staphylococcus aureus infections. PLoS One. 2009, 4 (5): e544610.1371/journal.pone.0005446.PubMed CentralView ArticlePubMedGoogle Scholar
 Ramilo O, Allman W, Chung W, Mejias A, Ardura M, Glaser C, Wittkowski K, Piqueras B, Banchereau J, Palucka A: Gene expression patterns in blood leukocytes discriminate patients with acute infections. Blood. 2007, 109 (5): 20662077. 10.1182/blood200602002477.PubMed CentralView ArticlePubMedGoogle Scholar
 Chaussabel D, Quinn C, Shen J, Patel P, Glaser C, Baldwin N, Stichweh D, Blankenship D, Li L, Munagala I: A modular analysis framework for blood genomics studies: application to systemic lupus erythematosus. Immunity. 2008, 29: 150164. 10.1016/j.immuni.2008.05.012.PubMed CentralView ArticlePubMedGoogle Scholar
 Xiao P, Chen Y, Jiang H, Liu Y, Pan F, Yang T, Tang Z, Larsen J, Lappe J, Recker R: In Vivo GenomeWide Expression Study on Human Circulating B Cells Suggests a Novel ESR1 and MAPK3 Network for Postmenopausal Osteoporosis. Journal of Bone and Mineral Research. 2008, 23 (5): 644654. 10.1359/jbmr.080105.PubMed CentralView ArticlePubMedGoogle Scholar
 Soetaert K, Petzoldt T: Inverse modelling, sensitivity and monte carlo analysis in R using package FME. Journal of Statistical Software. 2010, 33 (3): 128.View ArticleGoogle Scholar
Copyright
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.