Keywords: Bone remodeling, Mathematical model, Basic Multicellular Unit, Allosteric Control, Michaelis-Menten

To help to understand the modelling process that occurs when a scaffold is implanted it is vital to understand the rather complex bone remodelling process prevalent in native bone.

We have formulated a mathematical model that predicts osteoactivity both in scaffolds, as well as in bone in vivo and could set a basis for the more detailed allosteric models. The model is extended towards a bio-cybernetic vision of basic multicellular unit (BMU) action, when some of the regulation loops have been modified to reflect the allosteric control mechanisms, developed by Michaels-Menten, Hill, Koshland-Nemethy-Filmer, Monod-Wyman-Changeux.

By implementation of this approach a four-dimensional system was obtained that shows steady cyclic behaviour using a wide range of constants with clear biological meaning. We have observed that a local steady-state appears as a limiting cycle in multidimensional phase space and this is discussed in this paper. Physiological interpretation of this limiting four-dimension cycle possibly related to a conservative-like value has been proposed. Analysis and simulation of the model has shown an analogy between this conservative value, as a kind of substrateenergy regenerative potential of the bone remodelling system with a molecular nature, and to the classical physical value -energy. This dynamic recovery potential is directed against both mechanical and biomechanical damage to the bone. Furthermore, the current model has credibility when compared to the normal bone remodelling process. In the framework of widely recognised Hill mechanisms of allosteric regulation the cyclic attractor, described formerly for a pure cellular model, prevails for different forms of feedback control. This result indicates the viability of the proposed existence of a conservative value (analogous to energy) that characterises the recovery potential of the bone remodelling cycle. Linear stability analysis has been performed in order to determine the robustness of the basic state, however, additional work is required to study a wider range of constants.

A C C E P T E D M A N U S C R I P T

Introduction

The understanding of the mechanism of bone tissue regeneration, including metabolic and cellular regulation, has a clear practical application in tissue engineering and is therefore of great interest to both clinicians and researchers. Moreover, bone is of particular interest as only two types of active cells participate in bone regeneration (commonly referred to as remodelling or turnover). It is widely thought that there are osteoclasts and osteoblasts which are active "participants" of the widely accepted concept of the basic multicellular unit (BMU) [START_REF] Compston | Bone Marrow and Bone[END_REF], which balances calcium homeostasis with skeletal modeling and repair. It is the simplest model of tissue regeneration in animals, so the study of it could generate approaches to understanding more sophisticated metabolic and cellular cycles which function in the body. Taking another perspective, because bone is the hard tissue in the body, which maintains its shape, there is a high risk of mechanical injuries. This adds an important social and medical dimension to the investigation of the bone remodelling physiology. Indeed, the bone implant market in the USA is worth many billions of dollars. Metal implants, used traditionally to replace bone, frequently suffer from adverse effects and nowadays one can see a trend to develop new approaches and technologies for bone implant manufacturing.

Special emphasis should be paid to the potential to apply rapid manufacturing technology to generate bio-resorbable scaffolds for tissue engineering, which has been recognised for several years [START_REF] Braddock | Born again Bone, Tissue Engineering for Bone Repair[END_REF][START_REF] Xiong | Fabrication of porous poly(L-lactic acid) scaffolds for bone tissue engineering via precise extrusion[END_REF][START_REF] Vozzi | Fabrication of PLGA scaffolds using soft lithography and microsyringe deposition[END_REF]. However, despite some promising results [START_REF] Taboas | Indirect solid free form fabrication of local and global porous, biomimetric and composite 3D polymer-ceramic scaffolds[END_REF][START_REF] Yan | Layered manufacturing of tissue engineering scaffolds via multi-nozzle deposition[END_REF][START_REF] Koegler | Osteoblast response to PLGA tissue engineering scaffolds with PEO modified chemistries and demonstration of patterned cell response[END_REF] the full potential of this approach can only be fully realised when the bone remodeling process is well understood.

Moreover, if this knowledge is incorporated into a mathematical model of the basic biological (biochemical and biomechanical) remodelling processes this could then provide an important tool to enable the design of scaffolds to be optimised. The current understanding of bone remodelling processes is based on certain assumptions, one of which is mentioned above concept of the BMU [START_REF] Compston | Bone Marrow and Bone[END_REF], in that the emerging activity is controlled by a number of feedback

A C C E P T E D M A N U S C R I P T

loops, including genetic, physiological and immune, which function at the tissue, cellular and molecular levels. On the last level, the participation of many molecular messengers is very difficult to investigate in vivo. Even the general animal semantic model of such processes is not yet completed and is still in the state of re-verification and continual refinement. One explanation for this could be the experimental difficulties in the measurement of the activity very tiny molecular messengers, with short life-spans, which are subjected to considerable binding constants to receptors (even when dead). Moreover, the difficulties of conducting biochemical experiments in vivo are magnified when working within the hard tissue environment. This obstacle forces researchers to develop cellular level models combining them with some participation of molecular messengers involved in the regulation loop.

Recently in [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF], an interesting mathematical study was performed which indicates that the paracrine and autocrine relations are very important parts of the cellular model of the bone remodelling cycle. However, different paracrine mechanisms could be involved and the approach of present work differs from study [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF] in an important molecular respect that we tried to study the possibility to change the focus by giving cell relationships based on the formal allosteric form, starting from Michaelis-Menten (M-M) one-site molecular control and extending the study to the well known Hill, Monod-Wyman-Changeux (MWC) and Koshland-Nemethy-Filmer (KNF) patterns. All these forms of control have the potential to represent allosteric regulation and this could be very interesting when studying the RANKL/OPG balance regulation, for example [START_REF] Lemaire | Modelling the interactions between osteoblast and osteoclasts activities in the bone remodeling[END_REF]. Classical Michaelis-Menten kinetics [START_REF] Michaelis | Die kinetik der Invertinwirkung[END_REF] is one the most working approximation of many models in different fields of biochemistry, microbiology and biotechnology, for example, in pharmacological models [11], chemostat models [START_REF] Lenas | Coexistence of three Competing Microbial populations in a chemostat with periodically varying dilution rate[END_REF], or batchkinetics models [START_REF] Tohyama | Modelling of the mixed culture and periodic control for PHB production[END_REF][START_REF] Srinivasan | Dynamic optimisation of batch processes I. Characterisation of the nominal solution[END_REF][START_REF] Smets | Optimal adaptive control of (bio)chemical reactors: past, present and future[END_REF]. A number of research publications discuss the Michaelis-Menten control approach applied to the enzyme network [START_REF] Heinrich | Efficiancy and design of simple metabolic systems[END_REF][START_REF] Hofmeyr | Taking enzyme kinetics out of control; putting control into regulation[END_REF][START_REF] Fell | Physiological control of metabolic flux: the requirement for multisite modulation[END_REF][START_REF] Elsner | Metabolic control analysis: Separable Matrixes and Interdependence of Control Coefficients[END_REF][START_REF] Ortega | Optimal Metabolic Control Design[END_REF][START_REF] Yildirim | Parameter estimation of nonlinear models in biochemistry: a comparative study on optimization methods[END_REF]. Recently, Michaelis-Menten kinetics has been used to describe the changing rates of cellular activity during bone resorption [START_REF] Martin | Sensitivity analysis of a novel mathematical model identifies factors determing bone resorption rates[END_REF]. At

A C C E P T E D M A N U S C R I P T

the same time, there are models discussed with respect to modelling of the molecular feedback control in ligand-receptor regulation and in ligand transport regulation [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF][START_REF] Maalmi | Ligand diffusion and receptor mediated internalization: Michaelis-Menten kinetics[END_REF][START_REF] Rattanakul | Modeling of bone formation and resorption mediated by parathyroid hormone: response to estrogen/PTH therapy[END_REF]. In our study Michaels-Menten control has been chosen, as a first-stage of the allosteric control extension of our cellular model [START_REF] Moroz | Phenomenological model of bone remodeling cycle containing osteocyte regulation loop[END_REF]. Based on this result we have adapted the model to employ other well known molecular control mechanisms (Hill [26,[START_REF]Technological and Medical Implications of Metabolic Control Analysis[END_REF], KNF [START_REF] Koshland | Comparison of experimental binding data and theoretical models in protein containing subunits[END_REF], MWC [START_REF] Monod | On the nature of allosteric transitions: a plausible model[END_REF]). In the framework of widely recognised models/mechanisms of allosteric regulation [START_REF] Michaelis | Die kinetik der Invertinwirkung[END_REF] the cyclic attractor, described formerly for a pure cellular model [START_REF] Moroz | Phenomenological model of bone remodeling cycle containing osteocyte regulation loop[END_REF], prevails for different forms of feedback control (M-M, Hill, KNF, MWC). This finding demonstrates the viability of the proposition of the existence of a conservative-like value (analogous to energy) that characterises the recovery potential of the bone remodelling process. This result indicates that the robust behaviour of the model is maintained from the simple cellular level to the molecular biochemical level of regulation.

The main objective of this study was to investigate the potential to find paracrine and autocrine parameters (following [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF]) in the form of allosteric regulation. Reformulation of the model in terms of allosteric control could generate the intermediate model from a cellular to a biochemical one.

Model development

Our model development is based on well known assumptions that steoclastic and osteoblastic formation activities are coordinated in framework of BMU and this coordination effectively balances calcium homeostasis with skeletal modeling and repair. The models presented in [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF][START_REF] Lemaire | Modelling the interactions between osteoblast and osteoclasts activities in the bone remodeling[END_REF][START_REF] Moroz | Phenomenological model of bone remodeling cycle containing osteocyte regulation loop[END_REF] are very useful since they provide a good basis for modeling of osteocell's spectrum behaviour. However, although our initial model [START_REF] Moroz | Phenomenological model of bone remodeling cycle containing osteocyte regulation loop[END_REF] and also the original model described by Komarova et al. [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF] predict many different modes of dynamic behaviour of the BMU in bone remodeling control, a number of limitations to this model were identified by the authors.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

6 Among these limitations is a need for improvement to the autocrine regulation loop function and at the same time the paracrine regulation loops employ quite a wide range of parameters (e.g. power ~ -0.5) that could be even beyond the biologically relevant range. Additionally, a number of publications indicate the importance of the level of osteocyte regulation [START_REF] Burr | Calculating the probability that microcracks initiate resorbtion spaces[END_REF][START_REF] Martin | Towards unifying theory of bone remodeling[END_REF][START_REF] You | Ultrastructure of the Osteocyte Process and its Pericellular Matrix[END_REF], the role of the osteocyte apoptosis as a part of the mechanotransduction control mechanism [START_REF] Noble | Bome microdamage and cell apoptosis[END_REF][START_REF] Taylor | The cell transducer in damage-stimulated bone remodeling: a theoretical investigation using fracture mechanics[END_REF][START_REF] Bonewald | Osteocyte biology: Its implications for osteoporousis[END_REF] and role of stress [START_REF] Nakamura | Suppression of osteoblastic phenotypes and modulation of pro-and anti-apoptic features in normal human osteoblastic cells under a vector-averaged gravity condition[END_REF]. The above issues have driven us to modify the initial model [START_REF] Moroz | Phenomenological model of bone remodeling cycle containing osteocyte regulation loop[END_REF].

Firstly, osteocyte's apoptosis in the bone remodeling regulation loop has been considered and secondly the autocrine and paracrine control has been enhanced to make it more biologically relevant. The autocrine and paracrine feedback function were chosen not in potential form but more akin to the ligand-receptor response/binding function (Hill, Michaelis-Menten, KNF or MWC). Such functions have the allosteric, competitive inhibition and other control degrees of freedom with a clear biochemical sense (rather than fractal values which are purely theoretical). We refined the regulation loops that control the activity of the BMU and attempted to introduce the cybernetic point of view, such that the control should be minimised from both the (catabolic) energetic point and metabolic point of view. For example, a reason for this could be the limitation of the transport into the bone of the energetic substrates such as ATP and oxygen, as well others substrates. Changes from the physiologically normal bone turnover rate could destabilise the metabolic optimality not only on the local (bone tissue) level but could also create a supply problem for the body as a whole. Consequently the total number of molecular messengers of the bone remodelling process should not exceed a particular limit.

There are some considerations that in order to produce a robust bone remodeling process, the regulation needs robustness at all levels of regulation hierarchy, and indeed cellular level. The molecular biochemical regulation loops at the tissue level are only organised by cell interactions via common compartments (body or tissue media) such as marrow-media, lacunocanalicular microcirculatory system of periosteocytic fluid in the case of bone. The

A C C E P T E D M A N U S C R I P T

active interaction units are the osteocells. Phenomenological cellular models must reflect the tissue infrastructure of regulation and function in the whole body, as well as the very robust biochemical pathways. In the case of a multicellular organism and tissue, an evolutionary process is taking place. Probably, from the evolutionary point of view, the cellular level is even more important, because the multicellular body evolved from cellular colonies with initially poor communication. Taking into account this point the authors have attempted to develop and analyse the possibility of a cellular model and robustness at this level. The resulting cell-level control scheme based on the introduction of the osteocytes (OCt) control loop could be presented as in Fig. 1, where OCt apoptosis initiates the osteoclasts (OCl) maturation from osteoclasts precursors (solid arrows). Alternatively, the terms paracrine and autocrine are just macroscopic formalisations of the action of microscopic local factors, resulting in a form of feedback control of the dynamic system model. The number of reported local factors is quite wide and shows the complexity of regulation at the micro-level.

It is known that OPG-RANK-RANKL pathway is the major pathway involved in bone remodelling control, see for example, [START_REF] Hofbauer | The roles of osteoprotegerin and osteoprotegerin ligand in the paracrine regulation of bone resorption[END_REF][START_REF] Kong | Molecular control of bone remodelling and osteoporosis[END_REF][START_REF] Theoleyre | Cellular activity and signaling induced by osteoprotegerin in osteoclasts: involvement of receptor activator of nuclear factor nB ligand and MAPK[END_REF]. The autocrine effect occurs due to pre-OCl expression of RANK which is targeted by RANKL. It is interesting that RANKL can exist as a soluble protein about 31kDa so that OBl/Stromal cells can induce OCl formation in the absence of direct cell-to-cell contact, [START_REF] Horowitz | Control of osteoclastogenesis and bone resorption by members of the TNF family of receptors and ligands[END_REF]. This forms another paracrine/autoctine degree of freedom in BMU regulation. Effectively, the RANK-RANKL balance is regulated by osteoprotegerin. The first incorporation of this pathway in a mathematical model was described in [START_REF] Lemaire | Modelling the interactions between osteoblast and osteoclasts activities in the bone remodeling[END_REF]. However, as it follows from a number of studies, many others hormones, like PTH, growth factors, cytokines, vitamins and ions are involved in autocrine/paracrine regulation of bone resorption and formation. There are many indications that the number of factors (BMP-3, BMP-7, IGF I, IGF II, TGF-β3, FGF-2, VEGF) expressed by osteocytes are involved in autocrine/paracrine regulation of OBl and OCl, Heino [START_REF] Heino | Conditioned medium from osteocytes stimulates the proliferation of bone marrow mesenchymal stem cells and their differentiation into osteoblasts[END_REF] reported on the MLO-Y4 OCt

A C C E P T E D M A N U S C R I P T

factor from conditioned media that supply the growth factor with stimulatory effect on OBl.

This suggests that OCt plays a fine-tuning role in bone remodelling. Bakker [START_REF] Bakker | The production of nitric oxide and prostaglandin E2 by primary bone cells is shear stress dependent[END_REF] who studied the possible role of nitric oxide in bone remodelling, suggests that OCt apoptosis attracts OCl thereby activating remodelling. Westendorf [START_REF] Westendorf | Wnt signaling in osteoblasts and bone diseases[END_REF] stressed the role of secreted 39-46 kDa cystein rich glycoproteins (Wnts) and their role in signalling in OBl. Moseley [START_REF] Moseley | Interleukin-17 family and IL-17 receptors[END_REF] discussed the role of Interleukin-17 family (IL-17) cytokines, secreted by T-cells and their role in cancer metastasis to bone and regulatory effects in OCl precursor maturation. Additionally, bone cells also express a wide range of the neurotransmitter receptors as glutamate, γ-aminobutyric acid, purines and pyrimidines. The overall complexity of the regulation pattern of these local and homeostatic factors participating in differentiation, maturation and osteocells activity is illustrated in Fig. 1.

Regardless of the fact that the effect of the majority of these factors is not direct and is mediated sometimes by a long sequence of other molecular intermediates/stages, the effective action could follow basic regulative forms, like those of Michaelis-Menten or Hill. For example, in the model incorporation of the OPG-RANK-RANKL pathway, described in [START_REF] Lemaire | Modelling the interactions between osteoblast and osteoclasts activities in the bone remodeling[END_REF],

PTH involvement is essentially nonlinear and rather Michaelis-Menten in nature. Furthermore, some research [START_REF] Rattanakul | Modeling of bone formation and resorption mediated by parathyroid hormone: response to estrogen/PTH therapy[END_REF] described an effective model based on the premise that the inhibitory effect on TGFb1 (TGFb1 -induced production of OPG by marrow osteoblasts stromal cells) reducing effectively RANKL accordingly to Michaelis-Menten kinetics, when the effective M-M constant for OCl activity was introduced.

The strategy adopted to further develop the model was to compromise between the level of microscopic interacting molecular factors and the macroscopic form of feedback function in the phenomenological model of regulation. Taking into account that in order to produce a robust bone remodelling process, the regulation needs robustness at all levels of control, and indeed at cellular and molecular levels. The final loop of regulation links together

A C C E P T E D M A N U S C R I P T

(in some way) the participating cells and their precursors, as well as the bone material and even integrates body homeostatic systems like ion balance or the immune system.

The full set of differential equations for the general dynamic model, that includes different types of feedback mechanisms into our cellular regulation model [START_REF] Moroz | Phenomenological model of bone remodeling cycle containing osteocyte regulation loop[END_REF] could be rewritten as

) ( ) ( ) , ( ) ( ) , ( ) , ( ) , , ( ) , ( 1 1 2 2 4 3 4 3 4 3 1 1 2 2 3 2 1 3 3 1 1 s B s B OCt OCt OBl FB OCl FB x x f x x f x x s f x f x x x f K x f x x x x x f K x f x x - - - = - = - = - = - + - + -
x is the relative population density of osteoclasts (OCl), 2

x is the relative population density of osteoblasts (OBl), 3

x is the relative population density of osteocytes (OCt), 4

x is the relative change in bone mass (B), s is the level of mechanical stress, and s x 1 , s x 2 -are the steady-state values of 1

x , 2

x variables, K 1 and K 3 -are the receptor dissociation constants.

Functions i i i i i FB M M x Km x Km x f + = - ) , ( , n i n i n i i i FB Hill x K x K x f + = ) , ( , 2 2 2 2 2 1 ) , ( i i i i i i i i i i FB KNF aKd x Kd x aKd x zKd x K x f + + + = , i i i i i i i i i i i i i i i i FB MWC K K x K x K K K x x K K x x K x f 2 ) 1 1 ( ) 1 1 ( 3 2 ) 1 1 ( 1 ) 1 1 ( ) , ( 2 2
+ + + + + + = ,
describe different types of non-linear feedback loops of regulation (where K i , Km i , a, z, K d , K1 i , K2 i , K3 i are constants), Michaelis-Menten one-site molecular control [START_REF] Michaelis | Die kinetik der Invertinwirkung[END_REF] and other well known Hill [START_REF] Hill | Introduction to Statistical Thermodynamics[END_REF], Monod-Wyman-Changeux [START_REF] Monod | On the nature of allosteric transitions: a plausible model[END_REF] and Koshland-Nemethy-Filmer [START_REF] Koshland | Comparison of experimental binding data and theoretical models in protein containing subunits[END_REF].

Functions - OCl f , - OBl f , - OCt f , - B
f describe negative regulation feedbacks due to apoptosis, death or transformation of cells. In our study we operated with relative popular densities or their changes of cells, normalised for a steady state to 1. We tried to normalise all constants (Michaelis-Menten, p50 concentration) to unity. All initial rate constants were chosen following Komarova [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF]. The feedback functions parameters were chosen so that p50 for all functions remained the same at relative concentration parameter equal to unit, Fig. 2.

An examination of a range of feedback control functions shows that it is possible to model diverse molecular mechanisms of regulation. It is well known that the molecular local factors act by a number of very specific molecular mechanisms, finally expressed in certain non-linearities in kinetic equations. Different allosteric ("other site") forms of factor-receptor regulation is an important molecular mechanism control of cell functions, and, particularly, cooperativity is the interesting degree of freedom in such a regulation, because it characterises the degree (sharp or gentle) with which the regulation reaches a threshold.

Our study leads us to suggest that it is more likely that the ideal cooperativity in the case of receptors exposed on osteocells, is limited to 2-3. This means that the number of binding centres on a receptor is therefore higher then 2. For example, TGF-b has two subunits, which bind to type-II and type-I receptors. Binding to the binding site on type-II receptor

A C C E P T E D M A N U S C R I P T

causes the receptor to recruit by binding to the second binding site on the type-I receptor [START_REF] Wingfield | Tumour necrosis factor is a compact trimer[END_REF][START_REF] Thompson | Beta A versus beta B: is it merely a matter of expression?[END_REF]. After phosphorylation of the type-I receptor it recruits and phosphorylates Smad2 or Smad3 in the long chain by targeting the TGF-B response element on DNA. In that way the entire process is cooperatively regulated. So Type-I and type-II TGF-β receptors are likely dimmers what could be associated with inhibitory factors and follow allosteric models of regulation having cooperative character. TNF-α factor is a trimer [START_REF] Wingfield | Tumour necrosis factor is a compact trimer[END_REF] and NF-kB regulatory receptor-activator protein it is likely that it has three binding sites exposed into extracellular medium. The existence of RANKL in a soluble form about 31kDa is an interesting explanation that OBl/Stromal cells could induce OCl formation in the absence of direct cell-to-cell contact [START_REF] Horowitz | Control of osteoclastogenesis and bone resorption by members of the TNF family of receptors and ligands[END_REF]. This is another argument for considering molecular models of binding in BMU regulation containing alosteric regulation.

However, a little quantitative experimental data is available for calculation cooperativity (the Hill constants) of ligand-receptor binding of any part of osteo-regulation BMU systems, even for the most studied OPG/RANK/RANKL regulatory pathway. We tried to find a quantitative confirmation for our suggestions by calculating the Hill constant for data adopted from available literature sources. We calculated an example for the data adopted from [START_REF] Hoare | Measurement of Agonist and Antagonist Ligand-Binding Parameters at the Human Parathyroid Hormone Type 1 Receptor: Evaluation of Receptor States and Modulation by Guanine Nucleotide1[END_REF] for hPTH, from [START_REF] Beaudreuil | Molecular characterization of two novel isoforms of the human calcitonin receptor[END_REF] for calcitonin, from [START_REF] Regmi | Suramin interacts with RANK and inhibits RANKL-induced osteoclast differentiation[END_REF] for RANKL. Results are shown in Fig. 2B. From data established for calcitonin [START_REF] Beaudreuil | Molecular characterization of two novel isoforms of the human calcitonin receptor[END_REF] we calculated Hill coperativity k=1.2 (R 2 =0.95). It is well known that for many biosystems cooperativity is greater than 2 (Hemoglobin, Erithrocrourin for example), but it is rather for self-regulated water-dissolved systems, than for ligand-receptor systems.

But though the number of binding sites (ideal cooperativity of a system) is greater than 2 or 3, real cooperativity (rather Hill cooperativity) could be much less. For example, for hemoglobin-oxygen homeostatic system the number of binding sites (ideal cooperativity) is 4, however, calculated Hill cooperativity is about 2.7 (according to Hill, 1962 [26]) under normal physiological conditions.

Taking this into account we limited in our study the number on binding sites on the receptor system to the range 1-3, which from our point of view overlaps the real range for many receptors.

Finally, the resulting system appears as a generalisation of the model [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF], molecular approaches from [START_REF] Lemaire | Modelling the interactions between osteoblast and osteoclasts activities in the bone remodeling[END_REF] and an extension of our own model [START_REF] Moroz | Phenomenological model of bone remodeling cycle containing osteocyte regulation loop[END_REF]. The dependence in these feedback functions is undertaken with the purpose of reducing the number of parameters based on the following assumptions:

+ OCl f -osteoclasts regeneration has autocrine properties and is initiated by osteocyte's apoptosis, following [START_REF] Noble | Bome microdamage and cell apoptosis[END_REF][START_REF] Taylor | The cell transducer in damage-stimulated bone remodeling: a theoretical investigation using fracture mechanics[END_REF][START_REF] Bonewald | Osteocyte biology: Its implications for osteoporousis[END_REF]; osteoclasts degradation function,

) 1 ( ) 1 ( 2 2 1 1 4 3 4 3 3 2 23 4 22 3 2 1 1 1 2 2 2 1 11 3 2 1 3 3 3 1 1 1 - + - - = - = - - - + = - - + = x k x k x sx x x x x x x K x x x x x x x K x x x n n n n n n & & & & α β β β α β β α (2)
where x 1 is the relative population density of osteoclasts; x 2 is the relative population density of osteoblasts; x 3 is the relative change in the population density of osteocytes; x 4 is the The first equation in system (2) describes auto-and paracrine OCl regulation of osteoclasts production (first term) and removal (last two terms) and follows Komarova et al.

A C C E P T E D M A N U S C
The difference from the Kamarova model [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF] in terms of OCl dynamics is that we put in a paracrine term (which is the control feedback loop from the osteocytes) based on suggestions that OCt play a key role in this regulation [START_REF] Noble | Bome microdamage and cell apoptosis[END_REF][START_REF] Taylor | The cell transducer in damage-stimulated bone remodeling: a theoretical investigation using fracture mechanics[END_REF][START_REF] Bonewald | Osteocyte biology: Its implications for osteoporousis[END_REF][START_REF] Nakamura | Suppression of osteoblastic phenotypes and modulation of pro-and anti-apoptic features in normal human osteoblastic cells under a vector-averaged gravity condition[END_REF]. The last negative term in the first equation describes OCl removal and reflects our proposition of the delayed paracrine control of OCl by OBl. In the second equation of the system, which concerns OBl, the first term describes the paracrine-like feedback control of OCl on osteoblasts whilst the other terms describe osteoblast transformation into osteocytes and their apoptosis/death. The major difference of system [START_REF] Braddock | Born again Bone, Tissue Engineering for Bone Repair[END_REF] from the Komarova model [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF] is the introduction of the third equation to the model, which describes osteocyte dynamics, where the first term is responsible for their transformation/differentiation from osteoblasts and second term describes the osteocyte apoptosis/death affected by the level of mechanical stress (s). The fourth equation of system [START_REF] Compston | Bone Marrow and Bone[END_REF] follows the last equation from the Komarova model [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF] but with the minor difference that the relative change was chosen as 1.0 instead of 100% as in the original model. The relative change was chosen as a variable to describe OCt population density, so all the variables describing OCl, OBl and OCt relative densities are dimensionless. This approach was adopted so that the cyclic modes of the model could be explored without specific values for the cell variables being required.

A C C E P T E D M A N U S C R I P T

Indeed, system (2) is scaled by using relative population densities. Time is scaled also, where τ=t/τ c is the dimensionless time of system (2) and τ c is the time scale of the order of hours. Taking τ c in this time scale allowed us to estimate the range order of rate constants.

Because we wished to explore the cyclic/periodic behaviour of the system in biologically relevant timescales (measured in hours, days or weeks) the smallest unit of time considered was 1 hour and this was used for τ c in our study. If constant k 1 from [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF], for example, is recalculated in cells -1 *hour -1 this changes the constant from 0.24%*cell -1 *day -1 to 0.0001 cell - 1 *hour -1 or to the dimensionless value 0.0001 which is actually in the middle range of this constant that we have studied: 10 -7 -0.1, Table 1, Fig. 10, diagram k1. The other rate constants, α 1 , α 2 , β 1 , β 2 considerably changed their sense and dimensions compared to the values employed by Komarova et al [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF]. Rate constant α 1 is no longer the paracrine feedback parameter for OCl-OBl but it becomes now the paracrine OCl-OCt feedback constant in our system. However, we retain its value in range of the original Komarova model. Taking the initial value for parameter α 1 as α 1 =0.002 hours -1 or ~3 cells*day -1 (see, for example, Fig. 4, Fig. 567.), which is a biologically realistic rate for OCl generation at a single remodelling site, the dimensionless range of values for this parameter is 10 -5 to 10 3 that covers 8 decimal orders.

A similar analysis and recalculation procedure was undertaken for all rate constants employed in the model and their dimensionless range is shown in Table 1, and in Fig. 10 when robustness has been studied.

Thus, regarding the variables used we believe that the relative population densities are dimensionless but the timescale of our resulting data is related to the order of hours. On the other hand, when the cycle period is scaled in certain time-units, the cycles in all phase planes are actually in the same time-units when the population densities are normalised to dimensionless units. Therefore rate constants could change values that would render the model less biologically realistic. However, the underlying behaviour of the model (within the range of

A C C E P T E D M A N U S C R I P T

constants used) is unaffected. We also suggested that the new rate constant α 3 that describes the rate of OBl transformation to OCt should be in the same range as OBl first order rate constant α 2 or even higher since the characteristic time of the bone regeneration could be shorter than the characteristic time of OBl production. In our study the α 3 rate constant varies for 6 decimal orders between 0.00001 and 10.0, what effectively covers the entire range of first order kinetic characteristic periods from 1 minute to several years. The characteristical biological cellular times are obviously in the range of hours to weeks or months. This range of values is based on the osteocells production rate suggested by Komarova -effectively the OBl production constant in [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF]. This value has been derived from the first order kinetics data from experimental histomorphometric data [START_REF] Parfitt | Osteonal and hemi-osteonal remodeling: the spatial and temporal framework for signal traffic in adult human bone[END_REF][START_REF] Kato | Establishment of an osteoid preosteocyte-like cell MLO-A5 that spontaneously mineralises in the culture[END_REF][START_REF] Vashishth | Determination of Bone Volume by Osteocyte Population[END_REF][START_REF] Skerry | Early strain-related changes in enzyme activity in osteocytes following bone loading in vivo[END_REF]. However, this is first order rate constant when actually the OBl equation is second order and sometimes even fractal [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF]. This illustrates that ultimately from first order kinetic data it is possible to just conclude the range of rate constants. Moreover, the experiments to study kinetic parameters in bone are really semi-invivo experiments because of the specific nature of hard tissue (bone), which prevents true invivo results from being determined. Gathering experimental data on the kinetic parameters of bone enzymes, receptors and other complex molecular structures and then to calculate and employ these rate constants in kinetic models is rather conjectural.

Numerical integration of the systems was done by a fourth-order Runge-Kutta subroutine rkfixed using MathCad®, MathCad 2000 Professional, MathSoft Inc., 1999.

Statistical calculations were performed using SAS v.9.1, procedure CANCORR.

Results

Our numerical study employs a particular range of functions (M-M; Hill n=2,3; KNF; MWC) in order to examine whether the model simulates the periodical modes of the bone turnover cycle behaviour. In a majority cases it produces the cyclic (periodic) behaviour within

A C C E P T E D M A N U S C R I P T

a wide range of constants, see Table 1 and Fig. 10. The timescale for model is not linked to the selected unit of time in a straightforward way because of the renormalization of the parameters of the model (Michaelis constants, dissociation constants in Hill case, p50 concentrations) to unity.

Evaluation of the initial constant values was based on experimental histomorphometric data [START_REF] Parfitt | Osteonal and hemi-osteonal remodeling: the spatial and temporal framework for signal traffic in adult human bone[END_REF][START_REF] Kato | Establishment of an osteoid preosteocyte-like cell MLO-A5 that spontaneously mineralises in the culture[END_REF][START_REF] Vashishth | Determination of Bone Volume by Osteocyte Population[END_REF]. One should note that the rate constants value, even for the linear systems similar to above Eq. ( 2), is not directly related to the value of the formation or resorption rates that could be measured in experimental conditions. Normalisation of the population densities/concentration parameters leads to redefinition of the values of rate constants in the model. In this case it is quite difficult to validate the model constants and timescale but some general results like the character of local steady state remains topologically the same.

Among the constants shown in Table 1, for particular constants shown directly in figures, we found that system (2) yielded a characteristic cyclic behaviour, Fig. 3A Comparing Fig. 3A, 3B one can see that the initial state system relaxes in the OCl -dimension to a plateau and trajectory "collapses" there to a steady for long period cycle. We couldn't find the relaxation explicitly similar to toroid-like as it was suggested for pure cellular model [START_REF] Moroz | Phenomenological model of bone remodeling cycle containing osteocyte regulation loop[END_REF] but the process clearly shows sometimes two different cyclic modes. For this particular case of limiting cycle the trajectories of the system can be interpreted biologically, taking into account the existence of a particular surface in four-dimension space (OCl, OBl, OCt and bone material) and following our findings from the cellular model [START_REF] Moroz | Phenomenological model of bone remodeling cycle containing osteocyte regulation loop[END_REF]. When constants affect the rate of decay of the long period model.

Applying the range of parameters that we studied to date, the trajectories have not demonstrated a tendency to behave as a strange attractor -contrary to the findings for the cellular model [START_REF] Moroz | Phenomenological model of bone remodeling cycle containing osteocyte regulation loop[END_REF]. One can possibly say that cyclic attractor prevails after the introduction of the Michaelis-Menten/Hill kinetics. The biological range of parameters was changed compared to the cellular model and in addition they were also modified by scaling and introduction of the Michaelis-Menten control and other constants like V max [START_REF]Technological and Medical Implications of Metabolic Control Analysis[END_REF].

From the simulations displayed in Fig. 3 for OCt and for the bone material effectively from -0.5 to 0.5. Thus, from the point of view of the variable ranges of attracting dynamical cycle the model is reliable, even when scaled to use relative densities of osteocells and bone material. Initially for recalculation for our model the rate constants were chosen in hours scale. After scaling, their value changed and apparently the timescale of our model changed accordingly. This means that we cannot state the exact unit for the period of the cycle and it thus more reasonable to limit our interpretation to dimensionless terms. However, we should stress that the cycle has two distinct periodic processes -one long period and a steady short period cyclic attractor. The short period cycle continues after the long time-period has decayed when β 11 and β 23 are quite high. That means that the system relatively quickly returns back to the steady state after it is disturbed.

The feedback function

i i i i i FB M M x Km x Km x f + = - ) , (
is not the only form that contains degrees of freedom for allosteric regulation. There is large number of possible logistic-like feedback control functions that model allosteric regulation, for example kinetics of ligand-binding in an oxygen transport system. Various other control mechanisms (also allosteric in nature) which operate at the molecular level also have an affect at the cellular level. The comparison of findings for different feedbacks indicated that the cyclic attractor prevails irrespective of the different feedback function, Fig. 3A,5B,7A,8. The OCl-OBl phase curve became the most asymmetically-shaped for all feedback type functions. The dynamical range also changes. The range of periods and decay ratios remains in the same range as for the model based on M-M and Hill feedback function, Fig. 3B, 5A, 7B. The shape and amplitude of the attractor changes slightly, this could be caused by different cooperativity in the regulation as shown in Fig. 2B. Logistic saturation-like functions are well known in many regulatory networks for example, in neural networks as the neuro-somatic threshold function [START_REF] Demongeot | Positive feedback circuits and memory[END_REF], and it is an interesting subclass for molecular networks, including cellular-and-molecular regulated network such as the BMU. In addition, their employment in BMU models leads to cyclic attractor states that could explain some periodic modes and the existence of energy-like conservative values on which the robustness of the networks could be based. Allosteric models have additional degree of freedom and very clear, interpretable values from the point of view of molecular control, compared to fractal spaces. Although the MWC and KNF allosteric mechanisms are probably too complicated to be involved in receptor-mediator interactions in the case of osteocells communication in BMU, we have been studying formal introduction of such a sort of feedback functions in our model. This introduction does not change the topology of the cyclic steady state of our model, Fig. 8.

A C C E P T E D M A N U S C R I P T

Analysis of eigenvalues of the linerised matrix of a dynamical system is the usual method to evaluate the character of stability, see for example [START_REF] Woods | Heteroclinic tangles and homoclinic snaking in the unfolding of a degenerate reversible Hamiltonian-Hopf bifurcation[END_REF]. Our parametric robustness studies were performed to determine the sensitivity that influence the stability of the cycle (basic equilibrium state) in two different ways: [START_REF] Compston | Bone Marrow and Bone[END_REF] with respect to a partial alteration of parameters and ( 2) when all parameters were simultaneously varied. In first case all parameters were spanned in a wide range (along up to 6-8 decimal orders of magnitude). Results are shown in Fig. 10 where we limited our study to the Hill case with n=2. For cases n=1 , where one can see that the character of the equilibrium changes with variations of the constant. However in this model the "geometrical" effect of these constants seems to be negligible.

Thus, generally we could conclude that within the framework of local changes, ~100 percent in logarithmic scale (~1000% in linear), the partial variations do not significantly change the character of the equilibrium, that it is cyclic attractor. However in the case of variation of all rate constants at the same time the situation could be different.

In the second case study [START_REF] Braddock | Born again Bone, Tissue Engineering for Bone Repair[END_REF] was when all parameters were simultaneously and randomly varied. Because the number of parameters is quite large it was difficult to study regularly all combinations in the 13 dimensional parameter space. Our investigations were performed by Monte-Carlo method using MathCad software. We randomly generated 15000 combinations of exponentially distributed values for a uniform population in logarithmic coordinates, the 13-dimensional hypercube of all parameters within the limits indicated in Table 1. The results are shown in Fig. 11 for M-M and Hill (n=2,3). The equilibrium characteristics of the points could be classified as in Fig. 11F. The results show that within the random combination of parameters the main steady state is one when 2 characteristic roots are equal to zero, and two other roots are complex. It proves that for the majority (80%) of parameter combinations discussed, the equilibrium is cyclic in nature, which means that the limiting cycle is very robust.

From all figures one can see that the majority of equilibrium points designed by dots, when all characteristic roots have negative real part, which means that the equilibrium point is stable. However, one can see from the graphs that another stable equilibrium is when the equilibrium point is characterised by two negative real roots and by the pair of imaginaryconjunctive eigenvalues of the linearised matrix designed as "o", see Fig. 10F. It is well-known that a Hopf bifurcation can occur when the positive equilibrium loses stability.

The steady solution, designated as a "dot", is loosing the stability as a result that the pair of complex-conjugated eigenvalues

ω ξ λ λ i ± = 2 1 ,
of linearised matrix of system [START_REF] Braddock | Born again Bone, Tissue Engineering for Bone Repair[END_REF] cross over into the right half plane such that ξ>0, Fig11F. On these figures all spectrums of possible equilibrium in our case are specified. The designation of the equilibrium points are shown in upper-right corner of every characteristic scheme. However additional studies are required to investigate the bifurcation structure of the model in a wide range of constants.

Local sensitivity analysis is very good tool for estimation of the effect of a single parameter variation on the output of a dynamical system. However, because of large number of parameters [START_REF] Tohyama | Modelling of the mixed culture and periodic control for PHB production[END_REF] makes it ineffective to study the influence of every parameter when all others are fixed. We suggest that the best way to study effects of changing parameters may be to evaluate

A C C E P T E D M A N U S C R I P T

the generalised effect of variation of vector of parameters on the vector (set) of generalised outcome parameters. For this evaluation we employed canonical correlation analysis.

Canonical correlation analysis is usually effective when studying the relationship of two different sets (linear combinations) of parameters. We applied it to study the influence of rate parameters on the character of the equilibrium point, treating real and imaginary part of roots as an opposite set of parameters. The results are summarized in Table 1. The variables Var1-Var3 are the most representative linear combination of the equilibrium parameters and their correlations with the rate parameters are show in the table. One can see from the table that with the range of studying parameters for all models the most influenced are following: s, β 11 , k 2 and α 2 . This pattern is very typical for all types of feedbacks -Michaelis-Menten and Hill n=2,3. It is interesting that the first three parameters are the sort of dumping parameters (describing the removal of cells) in the system, which probably effect the transition of system from steady state with all negative roots to a equilibrium state with a cyclic behaviour. Anyway these two equilibriums are the most populated within the equilibrium modes shown in Fig. 10F.

Discussion

Earlier we stressed that the mathematical models of bone remodelling [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF][START_REF] Lemaire | Modelling the interactions between osteoblast and osteoclasts activities in the bone remodeling[END_REF][START_REF] Moroz | Phenomenological model of bone remodeling cycle containing osteocyte regulation loop[END_REF] predict various modes of dynamic behaviour of the BMU. In this study, the role of different feedback control functions has been investigated in the framework of an extended dynamic system. In this model the osteocyte regulation at the cellular level is formally introduced and the autocrine and paracrine regulations are choosen in an allosteric form. Regarding the first development, there are a many of indications that the osteocytes play a vital role in signalling mechanical damage [START_REF] Compston | Bone Marrow and Bone[END_REF][START_REF] Skerry | Early strain-related changes in enzyme activity in osteocytes following bone loading in vivo[END_REF][START_REF] Lanyon | Using functional loading to influence bone mass and architecture[END_REF][START_REF] Tomkinson | The death of osteocytes by apoptosis in human bone is observed following estrogen withdrawal in human bone[END_REF][START_REF] Noble | Bone loading changes the number and distribution of apoptotic osteocytes in cortical bone[END_REF][START_REF] Noble | Identification of apoptotic changes in osteocytes in normal and pathological human bone[END_REF][START_REF] Burger | Mechanotransduction in bone-role of the lacunocanalicular network[END_REF]. At the same time this form of model is also potentially applicable to the number of allosteric regulation controls types, for example by the dissociation constants K 1 and K 2 .

A C C E P T E D M A N U S C R I P T

Presented in this article phenomenological approach is that the formulation and interpretation of paracrine and autocrine control is made in terms of allosteric regulation rather than fractal form, when it is formulated as regulation over the degree parameters [START_REF] Komarova | Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone modelling[END_REF]. Such approach has an additional benefit -allosteric control degrees of freedom, clearly interpretable from the point of view of molecular control. It can be seen from Fig. 4C, Fig. 6B, 7C that the variations in K m (M-M model) and K D (Hill model), Eq. ( 2), significantly affect the parameters of the attracting cycle, which are related to an increase of the energy-like value of the remodelling potential. However, we need to bear in mind that the remodelling cycle in bone is not a conservative-like system and a flow of substrate and energy resource is required to maintain it. But in a homeostatic sense it probably could provide shorter delay in bone tissue recovery effect, than in the case if it does not exist. In this case our idea could be formulated in terms of losses, like in a dual problem of optimal control, when by means of the employing of conclude that a cyclic process is the optimal from the perspective of regulation. We also could say that any increase (Fig. 4, 6, 7C, 9) of its amplitude and period could be relevant to some physiological situations. The increase of phase amplitude of oscillations and the frequency of oscillating systems in physics leads to an increase in the energy of the process. The biochemical nature of the existence of such a cyclic attractor for a bio-process means supply of resource (e.g. oxygen, ATP, other feeding substrate) to this cyclic process and a change in the amplitude (Fig. 4, 6, 7C, 9) and frequency of oscillations of this cycle could reflect a strengthening or weakening of the supply of this resource. In the case when the phase-speed of the cycle is too fast, for example a Paget's-disease-like physiological situation, overfeeding of the remodelling cycle occurs. When the energy-like remodelling potential is underfeeding, the risk of osteoporosis physiological situations increases.

From the perspective of optimal control the metabolic expenses for the organism to The models considered here show just one subset of BMU models, within one combination of parameters that could have the attracting limit cycle in the range of parameters with biological meaning. However, our robustness study revealed that in a quite wide range of parameters there is a local minimum for the system (1,1,0,0), which has a cyclic-like steady state. This local minimum is relevant to the minimum for the cellular model and the model also has two cyclic modes with big differences in frequencies. One of these two cyclic modes maintains a steady state, whilst the other mode decays in time.

Due to difficulties in determining the precise biological principles in play to explain the behaviour of the model there are clear implications if one regards the process from a classical control theory perspective. One could reasonably compare the remodelling process to a mechanical process with a steady state cyclic behaviour, which is then subjected to an input disturbance. This system then also oscillates at a lower frequency, which decays until the systems returns to the original steady state oscillation like delay-coupled oscillators [START_REF] Sen | Collective dynamics of delay coupled limit cycle oscillators[END_REF].

Although the rescaling in the model reveals the geometry, a problem with interpreting the results from the model is that after renormalization the variables are dimensionless and thus the rates and the allosteric constants have lost their clear biological relevance to the values derived from experiments. However, we should say that even in the case of experimental measurements in vivo, these constants couldn't be precisely determined due to difficulties in conducting experimental procedures with hard tissue.

A C C E P T E D M A N U S C R I P T

Further developments of this model include the derivation of the equations for scaffold material, which employ specific constants for sorption and resorption. In such a way it could be possible to model the integration of resorbable implant into the bone. Such equations could include the material parameters, scaffold design (porosity) parameters via fractality of the bone scaffold, surface modification parameters and cell enhancement parameters. Our results also indicate that in the framework of the generalised model (Eq.1) it is possible to find other steady states with a sound biological/biochemical interpretation. Furthermore, the bifurcation between these states could model the major mechanisms of control of BMU and bone remodelling.

Conclusions

In this paper a model of the bone turnover has been developed. It extends the number of molecular and cellular control loops which are more realistic from a biochemical and osteohomeostatical perspective, since the effect of local factors is considered to participate actively in BMU regulation. This implementation was based on the assumption that this form of control is a vital part of bone remodelling, as it allows the completion of the natural loop of control of BMU initiated from the bone marrow. Simulations of the model demonstrate four-dimensional oscillations when employing the range of constant rates that can be interpreted biologically.

The results presented here show that basic steady state has the form topologically equivalent to a cyclic attractor in four-dimensional phase-space. Although the system is not conservative the survival of this cycle in four-dimensional "osteoclast-osteoblast-osteocyte-bone material" space indicates formally that there could be a conservative-like value that could characterise the cyclic attractor energetically. Biologically, or rather physically, this value could be interpreted as a "substrate-energy potential" of the bone remodeling system for bone recovery.

This potential could be associated with the continuously operating BMU and provides a measure of the recovery potential of the BMU following mechanical and biochemical damage Menten, n, number of binding cites, is equal to 2.0. B, the Hill (logit) plot of some important in bone regulation molecular ligands: data adopted from [START_REF] Hoare | Measurement of Agonist and Antagonist Ligand-Binding Parameters at the Human Parathyroid Hormone Type 1 Receptor: Evaluation of Receptor States and Modulation by Guanine Nucleotide1[END_REF] for hPHT, from [START_REF] Beaudreuil | Molecular characterization of two novel isoforms of the human calcitonin receptor[END_REF] for calcitonin and for RANKL from [START_REF] Regmi | Suramin interacts with RANK and inhibits RANKL-induced osteoclast differentiation[END_REF].

  regulation by osteoblasts, osteocytes and osteoclasts themselves; + OBl f - osteoblast activation loop positive feedback loop related to bone material density, -OBl f - osteoblast's transformation to osteocytes, lining cells and their apoptosis; + OCt f -osteocytes differentiation from osteoblasts depends on the bone material generated, -OCt f -osteocytes apoptosis is dependent on the stress attitude and the bone density; and, finally, + B f -the bone mass formation is dependent on the osteoblasts concentration and the bone resorption function -B f depends on osteocytes level. In the case of Michaelis-Menten (n=1) or Hill (n>1) forms of feedback function and the linear suggestions about the cell removal, system (2) transforms to

  figures, we found that system (2) yielded a characteristic cyclic behaviour, Fig.3A, Fig.4 with

β 11 and β 23 are

 23 very low ~0.0001 the surface indicates two types of cyclic behaviour with two periods and the long period cycle decays with time. For constants β 11 =0.01, β 23 =0.01 one can see only exponential relaxation, Fig.3B. The change of the rate parameters in Eq. (2) does not change the topology of cyclic attractor (Fig. 4A, Fig.10show that the attractor topology is quite robust when β 1 and α 2 change). In time the long period mode disappears and the short period mode becomes steady and shows limit cycle properties. Changes in the rate parameters (α and β)

  , Fig.4B, Fig.5B, Fig.7A one can see that the relative population density of OCl and OBl change in the range 0-4 and 0-1 respectively. For the OCt's and bone material we have chosen, following [8], the relative changes. So, it follows from the graph (Fig.3, Fig.4B, Fig.5B, Fig.7A) that the range for these parameters is -0.2 to 0.2

(

  Michaelis-Menten) and n=3(Hill) the results are very similar. A star indicates the fixed value of the parameter which yields a steady state as shown in Fig.3. Qualitatively the pictures remain the same for a wide range of parameters, including parameters shown in Fig.6A.One can see from Fig.10, diagrams for α 1 (a1), β 1 (b1) and β 11 (b11), that in the framework (range) of about 7 decimal orders of these constants the character of equilibrium changes a little. Just higher values of rate constant β 11 changes the nature of equilibrium from cyclic attractor to saddle-focus, Fig.10, diagram b11.From the group of rate constants of the second equation of system (2), the constants α 2 , β 2 , β 22 and β 23 , just parameter β 23 , describing the osteoblast removal, significantly changes the nature of equilibrium, Fig.10, diagrams a2, b2, b22, b23, accordingly. Adjustment of β 23 results in the attracting cycle collapsing and becoming a saddle-node equilibrium point. We could conclude it is likely that the robustness of the system depends on the rate of the OBl apoptosis -biological form of the cell removal.Changes in the values of constants of third equation effects to a little degree the character of equilibrium, see Fig.10, diagrams for α 3 (a3) and s. The curves are asymmetrical and indicate equilibrium changes from cyclic attractor to the stable point when the attractor collapses within the 6 decimal orders of changes.A C C E P T E D M A N U S C R I P TACCEPTED MANUSCRIPT 20Along with the variation of k 1 and k 2 rate constants in the second equation of system(2), the nature of equilibrium remains the same in 1-2 decimal order, see Fig.10, diagrams for k 1 and k 2 . The diagrams show anti-symmetrical character of dependence and little changes in quality of equilibrium within about 3 decimal orders with further changes to saddle-point equilibrium.The K 1 and K 3 group of parameters is very important because these parameters define the molecular signal "transfer function" and the order of binding (M-M, Hill n=2, Hill n=3) determines the molecular gain of these transfer functions. One can see that with the change/variation of dissociation constants K 1 and K 3 the character of the equilibrium point changes a little. With the increase of K 1 or K 3 (which means inhibition of ligand-receptor binding), the cycle rather collapses to the state (designed by dots), which means a stable attractor (see diagram for K 1 ). Inhibition/increasing by K 3 shows no differences in the nature of equilibrium point in 2 decimal orders, Fig 10, diagram for K 3 . Fig.6 shows similar allosteric effect of K 1 and K 3 dissociation constants on the cycle. It corresponds to Fig.10, diagram for K 1

  the cycle the bone homeostasis minimises the substrate-energy losses for bone remodelling by balancing the metabolic cost of regulation against (shortening response time to mechanical/aging damage) and physiological function of the skeleton. Concerning the role of OCt apoptosis, one can see from Eq. (2) and Fig.2that our suggestion is that the feedback function which regulates OCl response to OCt apoptosis is S-shaped (akin with many control feedback functions) and could have Michelis-Menten, Hill, MWC or KFN form. Our findings indicate the survival of a cyclic attractor (Fig.4-7, 9) in the multi dimensional phase space over a particular range of constants. In physics, for example, in the case of the classic onedimensional harmonic oscillator, the conservative surface is described by a circle in twodimensional phase space of the coordinate and its derivative. This circle represents the conservative value, i.e. the mechanical energy, that reversibly transforms from potential form to kinetic form. In our multidimensional OCl-OBl-OCt-Bone case, the 4-dimentional cycle reflects the existence of a certain interchange, i.e. the transformation from the "metabolic andA C C E P T E D M A N U S C R I P Tkinetic form" of the BMU into "bone material potential" which is supported in the steady-state by the blood supply of the substrate-energy resources.The important attribute of the model is the asymmetry of the cycle, Fig.4-9. This asymmetry indicates the complexity of the phase relation in the BMU, osteocytes, bone mineral and organic components. From the physical point of view employed above, a basic bone remodelling steady state turnover exists and all regulations in the direction of increasing or decreasing this level could be considered in the natural range of adaptability of bone remodelling, within the physiological activities of the body. Studying the model leads us to

  support this cycle are balanced with the need to recover the skeleton function in the appropriate physiological time. It is clearly apparent that the regulation of the cycle should be maintained, by the diverse loops of molecular control, in a very precise way because the metabolic losses to support the remodelling cycle have to increase with any change of these parameters. It is evident that within the framework of this quite phenomenological model, the role of the diverse A C C E P T E D M A N U S C R I P T molecular factors in bone regulation, such as receptors and mediators, the state of the membrane, and hormonal or genetic system, are difficult to include and discuss. The roles of these or any other molecular messenger or substrate remain the subject of broad discussion in the biochemical literature, even for the generalised animal model, and so the development of a mathematical model, based on the molecular level of regulation in the bone, awaits more precise biochemical and biophysical data.

A

  C C E P T E D M A N U S C R I P Tto bone. Consequently, in the framework of very-well known mechanisms and their models of allosteric regulation the cyclic attractor, formerly described for pure cellular model[START_REF] Moroz | Phenomenological model of bone remodeling cycle containing osteocyte regulation loop[END_REF], prevails for different forms of allosteric (Michaelis-Menten, Hill, MWC, KNF) feedback control. Furthermore, it becomes a more physiological shape in the meaning of molecular regulative parameters. These results support the suggestions of the existence of a conservative energy-like value that characterises the recovery potential of the bone remodelling process. The model could provide the basis for explanation of Paget's-disease-like physiological situation of overfeeding of the bone remodelling cycle. A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T

Fig. 1 .Fig. 2

 12 Fig. 1. Simplified scheme of the cell interaction pattern in the extended BMU. Continuous arrows represent the control loops of regulation in BMU on the cellular level. Dashed arrows represent molecular control pattern. NO, nitric oxide, PG, prostaglandin; TGFβ, TGFβ1, tumor growth factors; RANKL, receptor activator of nuclear factor κB ligand; OPG, osteoprotegerin; IGF, insulin-like growth factor; PTH, parathyroid hormone; 1,25(OH) 2 D3, vitamin D; M-CSF, macrophage-colony stimulating factor.

Fig. 3 .BFig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .

 3456789 Fig. 3. Phase trajectories for the system (2). The graphical matrix of phase plots, first raw-

Fig. 10 .Fig

 10 Fig.10. Local equilibrium point robustness studies -spanning of the system (2) parameters.The range is shown below the graph. Vertical axis is the x2 variable of system[START_REF] Braddock | Born again Bone, Tissue Engineering for Bone Repair[END_REF]. By star the

Table 1 .

 1 The range of constants used in study and correlation between constant rates and

			Michaelis-Menten		Hill, n=2			Hill, n=3
	Coefficient	Range	Var1	Var2	Var3	Var1	Var2	Var3	Var1	Var2	Var3
	α 1	0.0001-0.1	0.0133	0.0046 -0.0078 0.0135	0 -0.0055 -0.0012	0.0012	0.0093
	β 1	1-100	0.1097	-0.0024 -0.0793 0.1186 -0.0328	0.1335 0.0672	-0.0192	0.2475
	β 11	0.0001-0.1	-0.0713	0.4808 0.0658 -0.0363 0.3639 -0.0217 -0.1798	-0.025	0.5082
	α 2	0.0001-0.1	-0.0078	0.0441 0.2346	0.004	0.153	0.1153 0.0646	-0.0192	0.1948
	β 2	0.01-0.5	0.0208	0.0085 -0.0387	0.02	-0.0202	0.0141 0.0126	0.0044	0.0085
	β 22	0.01-0.9	-0.2693	0.1436 -0.1017 -0.2679 0.0564 -0.0295 0.0491	0.5348 -0.0312
	β 23	0.001-0.1	0.0644	-0.1419 -0.1686 0.0236 -0.1907 -0.0485 -0.111	0.0065	0.057
	α 3	0.0004-0.05	0.1293	0.0027 -0.0979 0.1683 -0.0358	0.1595 0.0685	-0.0058	0.2481
	s	0.05-0.5	-0.6317	-0.185 0.0535 -0.6152	-0.085	0.084 -0.8906	0.0731 -0.0478
	k 1	0.00001-0.001	0.104	0.0208 0.0312	0.1122	0.0228	0.0454 0.0176	-0.0406	0.1423
	k 2	0.001-0.1	-0.1782	0.1883 -0.253 -0.2519	0.028	0.0257 0.0505	0.6263	0.0645
	K 1	0.01-10.0	-0.0362	0.0426 0.0643 -0.0292 0.0808	-0.032 0.0078	0.0021 -0.0129
	K 3	0.01-10.0	0.0105	-0.0138 -0.0132 -0.0207 0.0208	0.0018 -0.0463 -0.0616	-0.098
		canonical variables summarising and characterising the robustness of local
		equilibrium; Var1, Var2, Var3 -canonical variables.

& & & & (1)where 1