Angela Madeo

Tomasz Lekszycki

Francesco Dell'isola

A continuum model for the bio-mechanical interactions between living tissue and bio-resorbable graft after bone reconstructive surgery

Keywords: ODEs governing growth/resorption, Coupling between mechanics and biological stimuli, Artificial bio-resorbable material, Continuum mixture model, Load-induced replacement of artificial material with natural bone tissue, Numerical simulations

A continuum model for the bio-mechanical interactions between living tissue and bio-resorbable graft after bone reconstructive surgery

Introduction, list of biological phenomena we want to describe and state of the art

In this paper a model based on well-established methodologies used in continuum mechanics is adapted to the description of the complicated bio-mechanical phenomena occurring in bones reconstructed with grafts of artificial bio-resorbable material. Some of the most important among these phenomena are the external load-driven resorption and synthesis phenomena. Of course, while a living bone may be resorbed or synthesised depending on the biological signal present in the bone tissue, artificial material can eventually be only resorbed. This is one of the more complex peculiarities of considered bio-mechanical systems which we aim to encompass in our model.

Motivation from medical practice

The use of bone substitute materials is nowadays a common practice in orthopaedic, jaw and skull surgery: therefore the determination of the optimal properties of these materials poses important challenges (see e.g. [START_REF] Greenwald | Bone graft substitutes: facts, fictions & applications[END_REF], [START_REF] Kraus | Treatment Strategies in Thoracolumbar Vertebral Fractures: Are there Indications for Biomaterials?[END_REF], [START_REF] Mano | Biodegradable and injectable polymeric matrix composites for hard tissue replacement: state of the art and recent developments[END_REF], [START_REF] Jarcho | Calcium phosphate ceramics as hard tissue prosthetics[END_REF]). Many bio-compatible materials of different chemical composition, form, macro-and micro-structure, mechanical and biochemical characteristics are available: some can be regarded as more or less bio-degradable and/or bio-resorbable (see e.g. [START_REF] Park | A Bone Replaceable Artificial Bone Substitute: Morphological and Physiochemical Characterizations[END_REF], [START_REF] Hing | Bioceramic Bone Graft Substitutes: Influence of Porosity and Chemistry[END_REF], [START_REF] Schilling | Resorbability of bone substitute biomaterials by human osteoclasts[END_REF], [START_REF] Burchardt | The biology of bone graft repair[END_REF], [START_REF] Kim | Ceramic bioactivity and related biomimetic strategy[END_REF], [START_REF] Goto | Resorption of synthetic porous hydroxyapatite and replacement newly formed bone[END_REF]). It is evident first of all that the bone substitute materials used in surgery must be bio-compatible and must exhibit a sufficient mechanical strength: indeed the implanted graft serves after surgery, especially in case of big cavities in bone, as a supporting element. Another requirement which is being considered more and more important is their ability to integrate with a natural tissue. Actually, bio-mechanical characteristics and micro-structure of bone substitute materials must be conceived in order to promote their integration with natural bone tissue: in particular an appropriate porosity of the artificial material is necessary to enable living cells migration and their expected activities. A suitable size of pores is also required for the development of vascularization necessary for survival of bone cells.

It is clear that the best possible expected situation after graft implantation is the gradual replacement of artificial materials with a natural bone tissue, as this tissue has some capabilities of functional adaptation which are far to be shared, at the moment, by artificial materials. In the desired replacement process three different and simultaneous phenomena play a relevant role: bone and artificial material resorption and new bone tissue formation. They are not independent and their coordinated contribution in the remodelling process determine the actual topology changes and eventual variations of mechanical state of the reconstructed bone. The resorption of both bone tissue and artificial material and the formation of intracellular matrix are guaranteed by the action of so called "actor cells" which is triggered by a signal determined by the mechanical state of bone tissue. The interplay of resorption and formation processes determines the final properties of the reconstructed bone.

Some well-established biological facts on bone tissue resorption and synthesis and how we account for them in the presented model

In the present paper we pay particular attention to the last phase of bone healing in presence of bio-material when the bone tissue should replace the added artificial material. This process is possible because of the activities of cells specialized in mechanosensation, mechanotransduction and tissue resorption and synthesis, see e.g. [START_REF] Stephen | Bone mechanics handbook[END_REF], [START_REF] Sikavitsas | Review: Biomaterials and bone mechanotrasduction[END_REF], [START_REF] Nomuraa | Molecular events caused by mechanical stress in bone[END_REF], [START_REF] Ruimerman | A theoretical framework for strain-related trabecular bone maintenance and adaptation[END_REF]. In particular, we can start distinguishing two types of cells (which belong to the common class of so-called "actor cells"), namely the osteoblasts (specialized in new bone tissue formation) and the osteoclasts (which are able to resorb both natural bone and bio-resorbable material) which have different roles and which originate from two different kinds of precursor cells. These precursor cells share the capability of freely moving inside porous materials which have a suitably large interconnected porosity and they can then penetrate and deposit inside a pore in both the natural tissue and porous artificial material. We will assume that these two kinds of precursor cells are present everywhere in both the living bone tissue and artificial material under the condition that a suitable porosity is present. The precursor cells give rise to osteoblasts or to osteoclasts following the instructions of a signal generated by a third kind of cells called osteocytes: this signal is proportional to the deformation which the osteocytes can measure at a given point (for this reason they are also called "sensor cells"). In order to accomplish their task of measuring deformation energy at a given point, the sensor cells do not move and spend all their life in that point. Sensor cells originate from osteoblasts when these latter have accomplished their task of synthesizing new bone around them: when an osteoblast is completely surrounded by new natural bone tissue it changes its nature and becomes an osteocyte, i.e. a sensor cell which starts to measure deformation and to emit a signal proportional to its measured value. Some additional remarks on the transformation of precursor cells into an osteoclast or an osteoblast are needed at this point. In particular, if a precursor cell detects a very low signal it interprets this datum as a very low deformation which the cell associates to the fact that more bone than needed is present close to it. In this case the precursor cell deposits and becomes an osteoclast which starts resorbing the exceeding material surrounding it. Conversely, if a precursor cell detects a very high signal, it associates this measure of high deformation to the fact that there is no enough bone around it and this precursor cell deposits and becomes an osteoblast which starts synthesising new bone tissue. Surgery must assure the initial mechanical integrity and basic strength of reconstructed bone. After the first phases of healing, the natural process of remodelling tries to enable the formation of a bone tissue able to functional adaptation: a well-designed graft must favour this natural attitude. At this stage the graft is generally surrounded by a new tissue: however two different possible graft evolutions may occur subsequently. In a first one the graft and surrounding bone tissue micro-structure does not enable migration of actor cells into the artificial material far enough from the interface with tissue, or does not enable the vascular system growth inside the graft. In this case bone synthesis in the graft is possible only very close to the interface with the tissue and the majority of graft mass remains as an inclusion in the natural tissue. According to the second possible evolution, the precursor cells manage to penetrate and survive into the domain occupied by the artificial material. The osteocytes signals, coming from bone tissue, manage to activate osteoblasts inside the graft which start to synthesize bone tissue. The action of osteoblasts and osteoclasts may therefore fill the pores in a graft with a new tissue or may result into artificial material resorption and its eventual replacement by a new natural tissue. All of these processes are interrelated, depend on many factors and must be efficiently described. Indeed, bone tissue reformation in a region initially occupied by bio-resorbable material determines how good is the final quality of reconstructed bone tissues. In this study we focus our attention on this final stage of healing after graft implantation in bone as this process is crucial for the success or failure of surgery. We consider a graft of porous bio-resorbable material located in a region in contact with bone tissue. To render the mathematical formulation of the model tractable it is useful to introduce a mixture of two porous materials: the grafted bio-material and the natural bone tissue. The apparent mass density and porosity of the two components of the mixture is allowed to change in time and space, together with its mechanical state, associated activity levels and fraction of all considered cells. The mathematical description we formulate in this paper for the bone remodelling process in presence of bio-resorbable material is based on the following assumptions:

(i) Osteoblasts and osteoclasts fraction and activity are proportional to the intensity of the stimulus they gather from the osteocytes: they elaborate and integrate all signals -sent by surrounding sensor cells-which reach them.

(ii) Osteocytes produce a signal proportional to the density of strain energy in the region where they are located. This signal is instantaneously transmitted and its intensity decreases exponentially with the distance from its source.

(iii) The sensor cells (osteocytes) are initially located only in the living tissue region (not in the artificial material one). Moreover, even for further times, since the osteocytes take birth from the osteoblasts when these latter are completely surrounded by natural tissue, the presence of sensor cells at a given material point is possible only if there is a contemporaneous presence of natural bone tissue at the same point.

(iv) The bone tissue can be resorbed or synthesized, while the artificial material can be only resorbed.

(v) The number of the cells which are active in a given place depends on the "effective" porosity of host composite material. In particular, the minimal values of this number are associated with great porosity (porosity = 1) and with absence of porosity (porosity = 0). Indeed, when the host material is too porous then the precursor cells cannot deposit at the surface of the host material, simply because there is not enough surface which they can attach to. When porosity is too low then the precursor cells cannot reach the considered region of host material. We therefore assume that there exists an optimal value of porosity associated with the presence of the largest number of cells and the highest level of their activities.

Main features of proposed model and comparison with similar models presented in the literature

The model which we present in this paper has many contact points with some continuum models already present in the literature. However, to our knowledge, the way we propose for describing the coupling between mechanical and biological phenomena is novel for what concerns the following features: i) the form of source term in growth/resorption ODEs, as depending non-linearly on stimulus and porosity; ii) the non-local dependence of local stimulus on deformation energy field in the whole body. For what concerns the specific phenomena we want to describe in this paper (and as described carefully e.g. in [START_REF] Turner | Mechanobiology of the Skeleton[END_REF]), the mechanical and biological structure of bone tissue is very complicated. A complete and simultaneous description of all the phenomena occurring in its growth or resorption may lead to models too complex to be treatable. Hence, we attempt here their synthetic description by means of a model reaching a suitable compromise between capability of describing the phenomena we are interested in and manageability. In this paper we use e.g. Cauchy continuum theory to model the mechanical behavior of the natural-bone/artificial-material mixture. However, it is well established (see e.g. [START_REF] Yoo | Couple-Stress Moduli of a Trabecular Bone Idealized as a 3D Periodic Cellular Network[END_REF], [START_REF] Rubin | The TEM Characterization of the Lamellar Structure of Osteoporotic Human Trabecular Bone[END_REF], [START_REF] Jasiuk | Modeling of Bone at a Single Lamella Level[END_REF], [START_REF] Ostoja-Starzewski | Couple-stress moduli and characteristic length of a two-phase composite[END_REF] and [START_REF] Bouyge | A Micromechanically Based Couple-Stress Model of an Elastic Orthotropic Two-phase Composite[END_REF],) that the particular micro-structure of bone tissue which actually is a composite material, may oblige us, in order to suitably describe some important mechanical phenomena, to use second gradient continuum theories. Indeed, in the present paper we focus our attention on the description of the process of graft resorption in reconstructed bones, assuming that a continuum mixture theory is sufficient to model the mechanical behavior of the more complex composite material constituted by artificial material and bone tissue. After having introduced a simplified description of stress state of bone tissue and artificial material mixture we address the problem of the choice of suitable constitutive equations for stress. To this purpose, we follow some of the indication found e.g. in [START_REF] Pawlikowski | Studies on constitutive equation that models bone tissue[END_REF] and [START_REF] Hamed | Multiscale modeling of elastic properties of cortical bone[END_REF] and we introduce suitable elastic moduli for the considered mixture, refraining, at a first stage, from any multi-scale modelling. Indeed, a multi-scale modelling should generalize the considerations developed for a pure bone tissue e.g in. [START_REF] Inglis | Characterization of anisotropy in porous media by means of linear intercept measurements[END_REF] or [START_REF] Yoo | Modeling of trabecular bone as a couple stress continuum[END_REF] (where again the trabecular bone is regarded as a couple stress continuum). In this first investigation we also consider the bone tissue/artificial material mixture as linear elastic: the present model needs to be improved by following the results shown e.g. in [START_REF] Pietruszczak | A fabric-dependent fracture criterion for bone[END_REF] where a fabric-dependent fracture criterion for bone is introduced.

In this paper an interface between pure bone tissue and artificial material will be considered: we will be able to describe some scale and boundary layers at this interface. The considerations developed in [START_REF] Wang | Scale and Boundary Conditions Effects on the Apparent Elastic Moduli of Trabecular Bone Modeled as a Periodic Cellular Solid[END_REF] are of relevance in this context. We believe that the final application of our theoretical investigations is the design of improved surgical techniques. For this reason we accept here the spirit which guides the investigations of many "bio-engineers" as those presented e.g. in [START_REF] Pietruszczak | Modelling of bone-implant interaction[END_REF] (where the bone-implant interaction is studied), in [START_REF] Skoworodko | Preoperative planning and post-operative estimation of vertebroplasty using CT/CAD/CAE systems[END_REF] (where a pre-operative planning procedure is conceived and a post-operative estimation of the obtained results is presented), in [START_REF] Dietrich | On concurrent engineering and design of an intervertebral disc of lumbar spine[END_REF] (where concurrent engineering and design of prostheses and surgical techniques which account for the bone remodeling phenomenon are developed). Even if only purely qualitative results are presented in this paper, we are able to show that the model we present should be able to reduce the occurrence of artificial material inclusions in the completely healed reconstructed bones. We finally underline the importance of the phenomena of load transfer between bone tissue and artificial reinforcement materials used in bone reconstruction (see e.g. [START_REF] Piszczatowski | Load transfer between elastic hip implant and viscoelastic bone[END_REF]): the numerical simulations we have obtained show how the gradual artificial material replacement is influenced by the modalities with which mechanical strain energy is differently distributed between the pure bone tissue and the bone tissue/artificial material mixture.

Eulerian and Lagragian description of a solid mixture constituted by bone-tissue and bio-resorbable material

For a reference concerning the basic concepts of continuum mechanics, which we use in this and all the subsequent sections, we refer to [START_REF] Germain | Cours de Mécanique des Milieux Continus[END_REF], [START_REF] Salencon | Mécanique des milieux continus[END_REF] and [START_REF] Stephen | Bone mechanics handbook[END_REF]. The considered solid material is constituted by a mixture of bio-resorbable material and living bone-tissue. In order to specify the kinematical descriptors which we want to introduce to describe the evolution of considered bio-mechanical systems, it is suitable to consider the following (Eulerian!) fields, evaluated at position x ∈ E, where E is the three-dimensional euclidean space, and at instant t: • the (apparent) mass density of bone-tissue

ρ b : (x, t) ∈ E × [0, T ] → ρ b (x, t) ∈ R + ;
• the (apparent) mass density of the bio-resorbable bone-substitute material (in the sequel, we also use the shorter expressions "artificial material" or " bio-material")

ρ m : (x, t) ∈ E × [0, T ] → ρ m (x, t) ∈ R + ;
• the porosity of the bone tissue/bio-material mixture i.e. the Eulerian volume fraction which is not occupied by bone tissue or bio-material

ϕ(x, t) : (x, t) ∈ E × [0, T ] → ϕ(x, t) ∈ R + ;
• the (apparent) total mass density and the bone tissue mass fraction

ρ(x, t) := ρ b (x, t) + ρ m (x, t), ξ b (x, t) := ρ b (x, t)/ρ(x, t);
• the (adsorption-growth) stimulus, i.e. a scalar quantity which measures the activation signal collected at (x, t) by the actor cells

S : (x, t) ∈ E × [0, T ] → S(x, t) ∈ R + ;
• the volume strain-energy density i.e. the energy needed to deform the elementary volume of the mixture bio-material/bone tissue from a reference configuration to the current one

U : (x, t) ∈ E × [0, T ] → U (x, t) ∈ R + ;
• the fraction of the active sensor cells, i.e. the number of active osteocytes over the number of potentially available osteocytes in the REV centered at (x, t).

d : (x, t) ∈ E × [0, T ] → d(x, t) ∈ R + .
This scalar parameter d is assumed to be dimensionless (like porosity).

The kinematical descriptors we have just introduced seem to supply the simplest mathematical model able to describe an elementary bone adaptation/remodelling process in presence of artificially added material. We now introduce a reference configuration to describe deformation occurring in considered solid mixture as it is well-known that the strain deformation drives the transformation of precursor into actor cells. Sensor cells, after having measured the deformation energy, will convey to precursors cells a suitable stimulus for the formation of the needed actor cells in the whole bone tissue/bio-material mixture. We assume that it is possible to suitably introduce • a reference (or Lagrangian) "average" configuration C * ⊂ E; • an "average" placement field χ : C * × R → E which allows for the introduction of the current configuration of the bone/bio-material solid mixture:

C t := χ(C * , t);
• two referential mass density fields

ρ * b : C * × R → R + , ρ * m : C * × R → R + ; • the total Lagrangian mass density ρ * := ρ * b + ρ * m ;
• the deformation gradient , its determinant and the Green-Saint-Venant strain tensor ,

F := ∇χ, J := detF, 2G := F T F -I,
in such a way that the following equalities hold

J -1 ρ * b = ρ b • χ, J -1 ρ * m = ρ m • χ, J -1 U * = U • χ, J -1 ρ * = ρ • χ.
The introduced reference configuration is the domain where all fields describing the remodelling process under study will be defined. It is most suitable and must be preferred to the Eulerian configuration, which is a time-varying domain.

Once the placement field χ is introduced, it is straightforward to define, for any X ∈ C * and for any time t, the classical displacement field

u : (X, t) ∈ C * × [0, T ] → χ(X, t) -X.
Presently we refrain from any effort directed to introduce any "homogenization" procedure, producing, after suitable averages, an expression for deformation energy in terms of "microscopic" phenomena occurring at a smaller scale. Indeed, the field χ is introduced to describe the "average" placement and deformation of the solid mixture under consideration, and is assumed to "fully" determine the deformation energy and able to relate the introduced referential and spatial mass densities.

Bio-Mechanical Constitutive, Evolution and Balance Equations

We start specifying, by means of suitable constitutive equations, the interaction between sensor cells and actor cells (see [START_REF] Ruimerman | A theoretical framework for strain-related trabecular bone maintenance and adaptation[END_REF]) as mediated by stimulus intensity. Indeed, as discussed in the introduction, it is well-established (see e.g. [START_REF] Stephen | Bone mechanics handbook[END_REF] and [START_REF] Ruimerman | A theoretical framework for strain-related trabecular bone maintenance and adaptation[END_REF]) that sensor cells produce a signal which is sent to precursor cells surrounding them and which drives their behaviour. Another phenomenon which we want to account for in our model concerns the effect of biological actions on mechanical behaviour of the mixture of bone tissue and bio-material: these actions affects the (apparent) densities of both materials, the mixture porosity and the elastic coefficients of bone tissue. We model the listed cause/effect relationships by means of simple evolution equations, when dealing with resorption of both bio-material and bone tissue or when dealing with synthesis of bone tissue. On the other hand, we introduce suitable constitutive equations when dealing with the relationships among biological and mechanical states, signal and stimulus, or stimulus and remodelling. We feel necessary to comment here about the modelling assumptions which have driven the formulation of the presented model. We assume a standard continuum mechanical model in which the elasticity coefficients can slowly vary with respect to time. Therefore, the standard considerations about constitutive equations and their relationship with second principle of thermodynamics hold. We are aware, however, that some energy must be supplied to the system in order to be assured that the postulated evolution equations for Lagrangian mass densities are valid. We are also aware that a general thermomechanical formulation encompassing also biological phenomena is needed and has been attempted in the papers [START_REF] Lekszycki | Optimality conditions in modeling of bone adaptation phenomenon[END_REF], [START_REF] Lekszycki | Modeling of Bone Adaptation Based on an Optimal Response Hypothesis[END_REF], [START_REF] Lekszycki | Functional adaptation of bone as an optimal control problem[END_REF]. In the present paper, we content ourselves to base the evolution equations of Lagrangian mass densities on phenomenological considerations which are guided by experimental evidence.

Deformation energy, fraction of sensor cells, porosity and stimulus constitutive Equations

The Lagrangian deformation energy density is assumed to be a function of the Green Saint-Venant strain tensor G, of the apparent volume Lagrangian mass densities and eventually of the considered material particle X. Indeed, differently to what happens in standard continuum mechanics, the Lagrangian mass densities may evolve with time in the considered resorption/synthesis problem. These changes are assumed to be slow enough so that inertia effects are negligible. In other words, we assume the existence of a function U * such that

J -1 U * (G, ρ * b , ρ * m , X) = U • χ.
In the sequel we will consider the following particular case of this constitutive equation in order to perform our numerical simulation, suggested by the biological phenomena which we have in mind:

U * (G, ρ * b , ρ * m , X) = μG 2 + λ 2 (tr(G)) 2 , where μ = μ(ρ * b , ρ * m , X) and λ = λ(ρ * b , ρ * m , X). ( 1
)
The introduced elasticity parameters are usually taken of the form (see e.g [START_REF] Stephen | Bone mechanics handbook[END_REF] and [START_REF] Ruimerman | A theoretical framework for strain-related trabecular bone maintenance and adaptation[END_REF]):

μ = μ 0 +μ 1b ρ * b -ρ * 0b ρ max A b +μ 1m ρ * m -ρ * 0m ρ max Am ; λ = λ 0 +λ 1b ρ * b -ρ * 0b ρ max B b +λ 1m ρ * m -ρ * 0m ρ max Bm , or even μ = μ 2b ρ * b ρ max C b + μ 2m ρ * m ρ max Cm ; λ = λ 2b ρ * b ρ max D b + λ 2m ρ * m ρ max Dm , ( 2
)
where all the quantities which have not been introduced before are constants and, in particular, the constant ρ max stands for the mixture apparent density corresponding to vanishing porosity. Both these equations for elastic moduli as functions of densities are admissible and seem to be coherent with experimental evidence available (see e.g. in [START_REF] Rho | Young's modulus of trabecular and cortical bone material: Ultrasonic and microtensile measurements[END_REF]). We choose to use equations (2) in our numerical simulations. As for the fraction of sensor cells, we assume that the biological activity of a material particle of the bone/bio-material solid mixture depends on its content of osteocytes. Therefore, in order to describe at least some of the biological mechanisms occurring in the solid mixture of bone and bio-resorbable material we need to introduce a constitutive equation relating the bone-tissue and bio-resorbable material mass densities to the fraction of sensor cells and porosity. According to the assumptions (iii) and (v) we made in section 1.2, we set

d = d(ρ b , ϕ), ( 3
)
where the fraction of sensor cells d is assumed to be an increasing function (ceteris paribus) of the total mass of living tissue in the mixture bone/bio-material. We now want to consider a constitutive equation for porosity. In order to do so, we recall that the biological process leading to the activation of osteoclasts and osteoblasts involves the diffusion of their precursor cells through the bone-tissue. Physiologically, bone-tissue is porous and a fluid is filling its pores. This fluid transports the precursor cells. In order to be assured that the activation process of precursor cells, leading to the formation and action of osteoclasts and osteoblasts, is made possible one needs to be sure that the considered mixture is porous in order to allow for the needed precursor cell diffusion phenomena. Moreover, the actor cells need to have pore surface where to sit in order to start their action. Of course, in the following considerations we will assume that the pores inside the bone tissue/bio-material mixture have a "regular" shape: this means that the pores have not any fractal nature, their shape ratio is not too big, and so on. This is needed in order to be assured that porosity is a correct measure of the available pore surface where actor cells may deposit. Porosity is a kinematical parameter which therefore needs to be introduced if one wants to describe remodelling phenomena in living bone-tissues and also in bone/bio-material solid mixture. In this mixture, if porosity is too low (and pores have a "regular" shape!) then living cells will not be able to efficiently adsorb the resorbable material neither to form new bone tissue, as the available space will not allow for the action of a suitably large number of actor cells. If porosity, on the other hand, is too large then we do not have enough matter (i.e. bone tissue or bio-material) on which actor cells may deposit.

Also in this circumstance remodelling will not occur quickly enough or will not occur at all. To describe the mechanical phenomena which influence the porosity variation in considered solid mixture we introduce a constitutive equation for the Lagrangian porosity ϕ * = ϕ • χ of the type

ϕ * = φ(G, ρ * b , ρ * m , X), (4)
which links the "effective porosity" (effective for the deposit of actor cells) to the strain-tensor G and to apparent mass volume densities. Behind this constitutive assumption for the porosity there is an implicit "instantaneous local microscopic equilibrium" hypothesis. We assume that locally (i.e. inside any elementary volume of considered mixture and at a microscopic level) transient phenomena, which adjust the porosity when local deformation state and apparent mass densities vary, occur in characteristic times which are negligible with respect to the time scale at which happen the macroscopic evolutionary phenomena which we want to describe. Moreover, we will assume that the final porosity in a macroscopic material particle of the considered mixture depends only on the "average" micro-elasticity moduli and on the "total" mass of bone tissue and bio-material present inside it and is nearly independent on the "micro" distribution of mass inside it. This assumption, rather standard in mixture theory, seems well-grounded even if maybe a more detailed analysis about it is required, by using suitable and delicate homogenization arguments. This assumption allows for the formulation of a macroscopic theory formally self-consistent: all informations about the phenomena occurring at micro-level are assumed to be accounted for by the introduced constitutive equation. Let us now consider the constitutive assumptions we want to make on the stimulus emitted by the sensor cells, as already stated at point (ii) of section 1.2. We assume that, similarly to what occurs in physiological bone-tissue, such a stimulus is conveyed in every material particle of the bone/bio-material solid mixture by the surrounding osteocytes. Moreover, we assume that the signal stemming from a sensor cell in a given material particle of the bone/bio-material mixture is instantaneously transmitted (i.e. the transmission time scale is negligible when compared with the characteristic time of the remodelling phenomena) and that its intensity decays exponentially with the distance from the emitting sensor cell. We also assume that the signal sent by any sensor cell located in a given material particle is of the same intensity, so that the total intensity of the signal sent by considered material particle is proportional to d. Finally, we assume that all sensor cells elaborate a signal which is proportional to the deformation energy concentrated in the elementary volume in which they are located. In conclusion, we assume that

S(x, t) = Ct U (x 0 , t)d(x 0 , t) exp(-D x -x 0 )dx 0 -S 0 (x, t) ( 5 )
where the length D -1 represents the range of action of sensor cells and the value S 0 is an activation threshold (eventually depending on the position and time) the meaning of which will become clear later, when a resorption-synthesis law will be formulated. For detailed discussion of the biological evidence which is the basis of the just formulated constitutive equation we refer to [START_REF] Stephen | Bone mechanics handbook[END_REF] and [START_REF] Ruimerman | A theoretical framework for strain-related trabecular bone maintenance and adaptation[END_REF]. As for dimensions of the stimulus, it can be checked that, since d is dimensionless, then S turns to have the dimensions of a work [N m]. Expression (5) for the stimulus is easily rewritten in the Lagrangian configuration by means of a simple change of variables

S * (X, t) = C * U * (X 0 , t)d * (X 0 , t) exp(-D χ(X) -χ(X 0 ) )dX 0 -S * 0 (X, t), ( 6
)
where we set

J -1 U * = U • χ, d * = d • χ and S * 0 = S 0 • χ.

Evolution Equations for Remodelling and Balance Equation of Force

We are now ready to formulate a well-posed evolutionary problem for the introduced kinematical fields which, in our opinion, is able to catch the most important features of the remodelling process occurring in bone tissue after initial healing and in presence of bio-material as a bone tissue filler and support. Again, we believe that remodelling phenomena occur without any inertial effect. Therefore, the evolutionary equation for apparent densities will be simply first order ordinary differential equations with respect to time. Moreover, we will assume that no spatial gradient of apparent densities occur in the remodelling equation. Also this assumption maybe discussed: we simply assume that the process of remodelling depends on what is occurring outside the elementary volume element only by means of the stimulus and is not influenced by other non-local (e.g. purely mechanical) phenomena. In formulas, we assume that the total mass density rate depends on the local value of the stimulus

∂ρ * b ∂t = A b (ϕ * , S * ), ∂ρ * m ∂t = A m (ϕ * , S * ) ( 7 )
The remodelling-constitutive equations for A b and A m are assumed to account for different phenomena of biological and mechanical nature and for some geometrical circumstances: i) the surface available for resorption or synthesis "inside" the considered macroscopic material particle depends, as already discussed, on the "effective" porosity field; ii) the different properties of bone tissue and bio-material determine different resorption rates, due to the different effect of actor cells on these different material; iii) the rate of synthesis of bone tissue which is consequence of a positive stimulus is different from the resorption rate determined by a negative stimulus of the same amount. The previous considerations are consequence of the biological nature of the process of synthesis and resorption: in particular, osteoclasts adsorb at the same time (and presumably they are not able to be completely selective!) both the bone-tissue and the resorbable material while osteoblasts produce only bone tissue. In particular, we will assume that

A b (ϕ * , S * ) = A b (S * )H(ϕ * ), A m (ϕ * , S * ) = A m (S * )H(ϕ * ), (8)
where the functions A b and A m are piece-wise linear functions with different slopes for negative and positive values of stimulus (remark that, according to assumption (iv) of section 1.2, for the bio-material and positive values of the stimulus, A m must vanish, as there cannot be synthesis of bio-material) and the function H is designed in order to account for the influence of "effective" porosity on the biological activity of actor cells: when "effective" porosity is too large there is not enough material on which actor cells may deposit, when it is too small there is not enough free space in the pores to allow their mobility and deposit. We choose the shape of H in such a way that H = 0 for ϕ = 0 or ϕ = 1. In particular, the following forms for the functions H(ϕ), A b and A m have been chosen to perform numerical simulations:

H = k ϕ * (1 -ϕ * ), A b (S * ) = ⎧ ⎨ ⎩ s b S * for S * > 0 r b S * for S * < 0 , A m (S * ) = ⎧ ⎨ ⎩ 0 for S * > 0 r m S * for S * < 0 , ( 9
)
where k is a suitable constant parameters, s b and r b will be called synthesis rate and resorption rate for bone tissue respectively and r m will be called resorption rate for bio-material. The evolution equations [START_REF] Dietrich | On concurrent engineering and design of an intervertebral disc of lumbar spine[END_REF], because of constitutive equations ( 8) are non-linear and depend on the value of strain energy of the whole body by means of the integral operator which gives the stimulus. Also this integral operator is strongly non-linear. However, they do not depend on the gradients of Lagrangian mass densities. In a sense, these mass densities are internal variables, following the definitions given e.g. in [START_REF] Engelbrecht | Microstructure described by hierarchical internal variables[END_REF], [START_REF] Engelbrecht | Hierarchical internal variables reflecting microstructural properties: application to cardiac muscle contraction[END_REF], [START_REF] Maugin | Thermodynamics with internal variables. 1. General concepts[END_REF].

We remark that, because of the presence of the stimulus in the argument on RHS of equation ( 7), the presented theory considers internal variables the evolution of which depends non-locally on coupled fields. Moreover, we have some difficulties in calling "observable" only the displacement field of considered continuum: our Lagrangian mass densities can be observed by means of X rays, even if, following the standard definitions, they are internal. Concerning the mechanical equilibrium of considered body, we will assume that inertia effects are not of relevance when considering the applied external load and the time scale of remodelling process. Therefore, we will assume that the loads applied to the system are varying so that a quasi-static deformation process is occurring. Hence, the balance of forces in the bulk reads:

Div F • ∂U * ∂G = -b ext , ( 10
)
where Div is the usual divergence operator, the central dot indicates a simple contraction between two tensors of order greater or equal than one and b ext accounts for the external body forces. As for the boundary conditions, we assume traction type conditions on the external boundary ∂C * and continuity of forces (which also implies continuity of displacement) at any internal discontinuity surface Σ * inside C * . In formulas

F • ∂U * ∂G • N = f ext on ∂C * , F • ∂U * ∂G • N = 0 on Σ * , ( 11
)
where f ext stands for the external forces per unit area applied at the boundaries of C * , N is the unit normal vector to this boundary (or to the internal discontinuity Σ * ) and [|a|] := a +a -stands for the jump of the quantity a trough the discontinuity surface Σ * . No difficulty arises when one needs to consider some "kinematical" excitation applied to the body, and we will not detail this point.

A simple one-dimensional, linearised case

Let us now consider a reference frame (x 1 , x 2 , x 3 ). We assume that the displacement vector has only one non-vanishing components along the x 1 axis, i.e. u(x 1 , t) = (u(x 1 , t), 0, 0). With this assumption, it is easy to show that

F s = ⎛ ⎜ ⎜ ⎜ ⎝ (u + 1) 0 0 0 1 0 0 0 1 ⎞ ⎟ ⎟ ⎟ ⎠ , G = 1 2 ⎛ ⎜ ⎜ ⎜ ⎝ [(u + 1) 2 -1] 0 0 0 0 0 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎠ .
where we clearly denote by an apex the partial differentiation with respect to the space variable x 1 . Moreover, we also assume that the vectors of applied external bulk and surface forces are of the form b ext = (b ext , 0, 0) and f ext = (f ext , 0, 0) respectively and we remark that in the considered one-dimensional case the unit normal takes the simple form N = (1, 0, 0). It can be recovered that in the introduced one-dimensional case the equation of motion and boundary condition given in Eqs. ( 10) and ( 11) take the following form

(λ + 2μ) 2 (u + 1) (u + 1) 2 -1 = -b ext on C * , ( 12
) (λ + 2μ) 2 (u + 1) (u + 1) 2 -1 = f ext on ∂C * , (λ + 2μ) 2 (u + 1) (u + 1) 2 -1 = 0 on Σ * .
where, in virtue of the constitutive equations ( 2) for the Lamé parameters λ and μ, we explicitly note that these parameters are not constant in space and time but depend on the Lagrangian densities of natural bone and artificial bio-resorbable material. These equations are then linearised around an undeformed equilibrium configuration thus simplifying into

(λ + 2μ) 2 u = -b ext on C * , ( 13
) (λ + 2μ) 2 u = f ext on ∂C * , (λ + 2μ) 2 u = 0 on Σ * ,
where, with a slight abuse of notation, u indicates here the perturbation of displacement around the equilibrium configuration. These one-dimensional, linearised equations will be used in the numerical simulations which we propose in the following section. As it is well known, these equations can be equivalently rewritten in terms of Young's modulus E and Poisson's coefficient ν as

(E u ) = -B ext on C * , ( 14
)
E u = F ext on ∂C * , [|E u |] = 0 on Σ * ,
To obtain these equations we consider that the Poisson's coefficient is constant and we set

B ext = b ext (1 + ν)(1 -2ν)/(1 -ν) and F ext = f ext (1 + ν)(1 -2ν)/(1 -ν).
We do not repeat here equations [START_REF] Dietrich | On concurrent engineering and design of an intervertebral disc of lumbar spine[END_REF] as they involve internal variables (see e.g. [START_REF] Engelbrecht | Microstructure described by hierarchical internal variables[END_REF], [START_REF] Engelbrecht | Hierarchical internal variables reflecting microstructural properties: application to cardiac muscle contraction[END_REF], [START_REF] Maugin | Thermodynamics with internal variables. 1. General concepts[END_REF]), i.e. variables the evolution of which is governed by ODEs with respect to time, so that their 1D version has a form which is not simplified with respect to their 3D version. The evolution equations for apparent Lagrangian mass density have a mathematical structure similar to the one studied in [START_REF] Luongo | Perturbation Methods for Bifurcation Analysis from Multiple Nonresonant Complex Eigenvalues[END_REF], [START_REF] Luongo | Multiple Scale Analysis for Divergence-Hopf Bifurcation of Imperfect Symmetric Systems[END_REF], [START_REF] Luongo | Multiple Timescales Analysis for 1:2 and 1:3 Resonant Hopf Bifurcations[END_REF], [START_REF] Di Egidio | Linear and non-linear interactions between static and dynamic bifurcations of damped planar beams[END_REF].

One must remark that, unfortunately, the analysis which leads to a "bifurcation chart" of the type found in figure 4 of [START_REF] Luongo | Perturbation Methods for Bifurcation Analysis from Multiple Nonresonant Complex Eigenvalues[END_REF] cannot be easily repeated in the present instance. Indeed, in our problem, the evolution equations for internal variables depend on the stimulus, which is itself obtained by an integral operator from the solution of the PDEs representing mechanical balance of force. Therefore, the study of bifurcation in the present instance also must account for the possible variations of resorption and synthesis parameters, which will be seen in the next sections to dramatically influence the final form and structure of reconstructed bone. As a consequence, we are forced, to rely on the result of numerical simulations in order to study how introduced constitutive parameters influence the final structure of reconstructed bone.

Non-dimensional form of mechanical and biological equations

We now want to write all the equations in non-dimensional form. To do so, we start from mechanical equations and we introduce the non-dimensional variables x = x/L, ũ = u/u 0 , Ẽ = E/E b , where L is a characteristic length of the specimen, u 0 is a characteristic displacement and E b is a constant Young modulus which we chose to be that of compact bone. With these notations, the dimensionless form of equations ( 14) reads

Ẽ ũ = - B ext L 2 u 0 E b =: -Bext on C * , ( 15
) Ẽ ũ = F ext L u 0 E b =: F ext on ∂C * , Ẽ ũ = 0 on Σ * ,
where, with a slight abuse of notation, the apex now indicates the derivative with respect to the non-dimensional variable x. In the sequel we will let the Young modulus of the mixture vary with density according to

E = E b (ρ * b /ρ max ) β + E m (ρ * m /ρ max ) β ,
E b = (λ 2b + 2μ 2b )(1 + ν)(1 -2ν)/(1 -ν) and E m = (λ 2m + 2μ 2m )(1 + ν)(1 -2ν)/(1 -ν)
. As a consequence, the corresponding dimensionless Young modulus thus takes the form

Ẽ = ρ * b ρ max β + Ẽm ρ * m ρ max β , ( 16
)
where Ẽm := E m /E b . Moreover, in the considered one-dimensional, linearised case the deformation energy given in Eq. [START_REF] Bouyge | A Micromechanically Based Couple-Stress Model of an Elastic Orthotropic Two-phase Composite[END_REF] simplifies into

U * = 1 2 (λ + 2μ)(u ) 2 = 1 2 (1 -ν) E (1 + ν)(1 -2ν) (u ) 2 ,
or equivalently, in terms of the introduced dimensionless variables

U * = U 0 Ẽ (ũ ) 2 , with U 0 = 1 2 (1 -ν)E b (1 + ν)(1 -2ν) u 0 L 2
We now want to find a dimensionless form of the biological evolution equations [START_REF] Dietrich | On concurrent engineering and design of an intervertebral disc of lumbar spine[END_REF]. To do so, we first introduce the dimensionless Lagrangian stimulus S * as

S * = S * U 0 L 3 = C * Ẽ (ũ ) 2 d * (X 0 , t) exp -D X -X0 + u 0 L ũ(L X) -ũ(L X0 ) d X0 -S * 0 , ( 17
)
where Eq. ( 6) has been used to get the second equality and where we introduced the dimensionless quantities D = DL, X = X/L, X0 = X 0 /L, C * = [0, 1] × [0, 1] × [0, 1], ũ = u/L and S * 0 = S * 0 /(U 0 L 3 ). Note that, since in linear elasticity the quantity u 0 /L is often very small, the third term inside the norm can be neglected in most of practical applications. Dividing equations ( 7) by ρ max , introducing the dimensionless time t = t/t 0 and considering the constitutive assumptions [START_REF] Goto | Resorption of synthetic porous hydroxyapatite and replacement newly formed bone[END_REF], their dimensionless form reads

∂ ∂ t ρ * b ρ max = H(ϕ * ) Ãb ( S * ), ∂ ∂ t ρ * m ρ max = H(ϕ * ) Ãm ( S * ). (18)
In these equations we set

H = k ϕ * (1 -ϕ * ), Ãb ( S * ) = ⎧ ⎨ ⎩ sb S * for S * > 0 rb S * for S * < 0 , Ãm (S) = ⎧ ⎨ ⎩ 0 for S * > 0 rm S * for S * < 0 , (19)
where

sb = U 0 L 3 t 0 ρ max s b , rb = U 0 L 3 t 0 ρ max r b , rm = U 0 L 3 t 0 ρ max , r m . ( 20
)
Finally we choose, for the particular numerical simulations we present in this paper, the following simple forms of the constitutive equations ( 3) and ( 4) for the Lagrangian fraction of sensor cells and for the Lagrangian porosity respectively

d * = η ρ * b ρ max , ϕ * = 1 -θ (ρ * b + ρ * m ) ρ max with, 0 < η ≤ 1, 0 < θ ≤ 1 (21)

Some numerical examples showing the potentialities of the presented model

In order to show that the simple model presented in this paper is actually able to describe different possible mechanisms of bone remodelling we address the study of a simple 1-D problem which seems to be the simplest possible presenting some of the required complex features of formulated model. In order to simplify the study of the mechanical part of the problem, we consider a sample which behaves as a rod (composed by two materials) with variable (with respect to time and space) axial rigidity and porosity. However, time variations of rigidities and variations of loads occur in a time scale which is much bigger than the transient characteristic time leading to macroscopic mechanical equilibrium. Therefore inertia effects can be neglected in evolutionary problems considered here. The numerical code used to calculate the sequence of configurations of the considered bar-bone includes a simple closed form expression which links external mechanical applied loads to the corresponding equilibrium deformation states, once the profile of space varying rigidity is given. In the present treatment only mechanical boundary conditions are needed. Therefore the evolutionary problem which we study will simply involve the equations for time variations of apparent mass densities of bone tissue and bio-material. In these equations the previously listed constitutive equations, and the parameters appearing in them, will appear, which establish the required link between mechanical and biological phenomena.

The used numerical code was written in C Language and we implemented a discrete time step increment algorithm. The flow chart which has been considered for formulating the numerical code used in the following section can be reconstructed when considering that it is composed by the following steps: (i) Calculate the variations of the values of ρ * b /ρ max and ρ * m /ρ max at time t by means of Eqs. [START_REF] Jasiuk | Modeling of Bone at a Single Lamella Level[END_REF]. N.B.: at first time step ( when we assume that the history of applied external loads starts to be known) the values of ρ * b /ρ max and ρ * m /ρ max are assumed to be also known (e.g. piecewise constant initial distribution of natural bone tissue and artificial material densities). (ii) Given the time history of externally applied load, use equilibrium equations [START_REF] Inglis | Characterization of anisotropy in porous media by means of linear intercept measurements[END_REF] (loads are assumed to vary quasi-statically) to calculate the displacement at considered time. (iii) For considered time, calculate fraction of sensor cells and porosity according to Eqs. and [START_REF] Lekszycki | Optimality conditions in modeling of bone adaptation phenomenon[END_REF] respectively. (iv) Calculate the stimulus at considered time according to Eq. ( 17). (v) Solve, for subsequent time step, the evolution equations ( 18) in order to get the value of the variation of the two densities at this subsequent time. (vi) Continue with step (ii). A very delicate issue in the performed numerical simulations regards the choice of time step. Indeed the integro-differential equations we are integrating are strongly non-linear. Therefore too large time step may direct the calculated evolution towards configurations which actually are non-reachable in reality. This issue may become a serious one when dealing with coupled PDE for mechanical equilibrium and the equations for bone tissue and bio-resorbable material time evolution. In the 1-D case which we will study in the next sections, however, a careful choice of time step depending on the values of constitutive parameters used is simply resolving the problem. We consider the case of a 1D rod sample of length L with initially completely separated bone tissue and bio-material thus assuming that one part of the rod is initially occupied by pure bone tissue and the remaining part is occupied by pure bio-material. The rod is clamped at x = 0 (u(0) = 0) and is loaded by a dimensionless static load F ext at x = L. In the case of compact bone, characteristic values of Young modulus and Poisson's coefficients are E b = 17GP a and ν = 0.3, so that, considering a characteristic displacement u 0 = 10 -4 L, a unit dimensionless force F ext = 1 corresponds to a dimensional applied force f ext = 2.3 10 6 P a. In other words, if the cross section of the specimen is 1mm 2 and the length is L = 10mm then one needs to apply a force of 2.3 N in order to get characteristic displacements of the order of u 0 = 10 -3 mm. We are interested into the phenomena occurring in the neighbourhood of the interface separating the two considered zones (it has to be remarked that most likely some interesting boundary layer phenomena may occur at this interface and that second gradient continua maybe of use in describing them, but we will postpone these considerations to further investigations, see e.g. [START_REF]dell'Isola and Seppecher, Commentary about the paper "Hypertractions and hyperstresses convey the same mechanical information Continuum Mech[END_REF]).

The effect of applied loads on the time evolution of healing process in surgically reconstructed bones in presence of bio-material graft

The following pictures illustrates two different possible scenarios depending on the level of externally applied forces. It is well known that loading affects dramatically the process of formation, healing and remodelling of bones. The biological system which naturally has been developed by evolution is capable of an action of remodelling excited by mechanical loading which is often called functional adaptation. Our aim is to determine under which conditions the surgical reconstruction with graft does not interfere or at least does not impeach the action of natural functional adaptation. More precisely, we want to determine in which conditions and to what extent the action of functional adaptation will lead to a replacement of bio-adsorbable material with bone tissue or at least to the final constitution of a stable composite. It has to be explicitly remarked that the presence of a very small bio-material apparent mass density is not forbidden when only "biological" processes are considered. Indeed when a small mass density for bio-material is present, the bone tissue from the interface may always penetrate into the region occupied by bio-material eventually reconstructing partially the rod bone. The problems related to void formations in presence of small density of bio-material are of mechanical nature: a weak bio-material graft is not able to sustain the applied mechanical loads! Therefore a warning should be formulated at this point: in our model it is not included the description of mechanical damage and crack formation. The material is always assumed to be elastic and in the performed simulations even linear elastic. Further planned investigation will include the inclusion of damage description and evolution in our model. As already pointed out in Tab. 1, we consider a constant initial distribution of the densities of bones and bio-material as shown in Fig. 1. In particular, we initially have only natural bone tissue on the left /ρ max = 0 and ρ * m /ρ max = 0.5) ). We remark that, due to the fact that we are using a two-solid macroscopic mixture theory, at each point of the specimen (x axis) one needs to have both a material point of natural bone tissue and of artificial material, eventually with different densities. We start discussing the time evolution of mass densities distributions for low level of mechanical loading shown in Fig. 2. In order to correctly interpret all the following figures, it must be remarked that the black colour stands for zero (or very low) concentration of the considered constituent (natural bone tissue or artificial material), the red colour stands for the natural bone tissue density and the blue color stands for the bio-material density. The intensity of the color is proportional to the density of the considered constituent, e.g. a more brilliant red means that the concentration of the natural bone tissue is higher. Thin bars display separately the distributions of bone tissue (red) and bone-substitute material (blue), while the thick bar reflects the sum of these distributions, so rendering the density distribution for the composite bone-tissue/bone-substitute material. In the authors opinion, figures of the type of Fig. 2 are suggestive, but we also report the same result in Fig. 3 in order to fit other stylistic tastes. We can observe in this case that an insufficient loading leads to a rapid resorption of bio-material everywhere, while bone tissue is resorbed in some regions and in others is formed. The formation of bone tissue is not rapid enough to fill the region of bio-material which is being resorbed too quickly (especially far from the interface with bone tissue) and the process results into an increasing void formation. These phenomena can be explained if one thinks that low external loads are associated to low value of the corresponding deformation energy and so to a low stimulus which initially stems from the region occupied by natural bone tissue. The precursor cells which are located far from the interface do not receive any signal and interpret this datum as very low deformation so transforming themselves in osteoclasts which start a quick resorption of the bio-material. The second considered case is illustrated in Fig. 4. The attained equilibrium configuration is also reported in Fig. 5 The described process has different features. As before, no (or very low) signal can reach the regions of bio-material far from the interface in the earlier times so that resorption starts occurring in these regions. The resorption process, however, increases the porosity so that the stimulus (which is much higher than the stimulus in the previous case due to a bigger value of the deformation energy measured by sensor cells located in the bone tissue) can reach the precursor cells located in the bio-material. These precursor cells thus become osteoblasts and start synthesising new bone tissue. At the end of the process, we can observe that the new bone tissue has almost completely replaced the bio-material far from the interface, while an inclusion of bio-material is left closer to the interface. This simulation hence predicts the presence of bio-material inclusions surrounded by bone tissues in the equilibrium configuration attained by grafted bones. We can explain the residual presence of bio-material by remarking that the porosity is not, in these inclusions, large enough to allow for the action of both kind of actor cells. Such kind of analysis may also enable a future optimization of both bio-material and mechanical conditions necessary to assure the absence of inclusions or, more generally, the final production of efficient bone tissues.

The effect of applied load on the final distribution of bone tissue and bio-material in rod-bones

The set of numerical simulations which we show in this section were aimed to consider the effect of external applied loads on the final result of the remodelling process producing the reconstructed bone after surgery. In the following figure we show the final result of remodelling process with varying externally applied loads. The simulations show that for different values of applied external mechanical load, at least two different ranges of values are possible. For lower values of external loads the resorption phenomenon is prevailing in both considered regions.

However in the region previously occupied by the bio-material there are no sensor cells so the only signal possible is the one originating from the sensor cells located in the bone tissue. This results into a lower stimulus and hence in a faster resorption in the region occupied by bio-material. Therefore, because of the absence of a sufficient number of actor cells being activated in this region, a void formation occurs. As the signal from the bone tissue towards the artificial region is too low, then the bio-material is finally bio-resorbed. Also in the bone tissue region the net result of the activity of actor cells produce resorption and also this tissue is finally reduced in mass. For higher values of externally applied loads a different sequence of events arises. The signal from bone tissue is much stronger and it penetrates into the region occupied by bio-material. This results into the following different steps in bone remodelling. First the bio-material far from the bone tissue starts to be resorbed, Therefore porosity is available for actors cells. Simultaneously close to the interface (and in the region initially occupied by bio-material) bone formation occurs resulting into the formation of a composite. These phenomena enables the increase of stimulus in the far region where actor cells resorbed the bio-material. In the so formed porosity, precursor cells become osteoblasts and new bone tissue appears and finally becomes dominant. Unfortunately when this process comes to an equilibrium configuration, porosity in the central part of the bar-bone reaches the values for which no actor cells can access the composite material and all biological activities stop. This is due to the fact that no resorption of bio-material is any more possible. The presented numerical simulations show how it may happen that in a remodelling process some bio-material may be residual. These results, although only qualitative in nature, may result into an interesting indication for practical surgery: the reconstructed bone must be loaded, but the load cannot be too large, as otherwise the residual amount of bio-material maybe too large, while in presence of too low loads large void formation occurs in the regions initially occupied by bio-material.

Conclusions

In this paper we formulate a mathematical model suitable for describing bone tissue synthesis and resorption in presence of a bio-resorbable material, to be used in bone reconstruction. Porosity of bone tissue and bio-material mixture is accounted for, as it plays a relevant role in the biological activity of actor cells which results into the resorption of both bone tissue and bio-material and into bone tissue formation. Different resorption and synthesis rates are considered possible in the presented model, together with different biological signal propagation distances and stimulus excitation. In order to investigate and explain some of the possible mechanisms of interaction of bone tissue and bio-material a simple one-dimensional case has been formulated and studied. To see how the different external loads influence the process of interaction between bio-material and bone tissue, in some numerical examples, we considered these loads as applied to bar-type bones, the dimensions of which are few millimeters. The performed numerical tests are so promising that we believe they must be used as a starting point for more general further investigations. In particular, it is shown that, in the framework of the presented model, one can forecast in which ranges of applied external loads, characteristic signal distance, synthesis and resorption rates and other bio-mechanical parameters, void formation occurs in the bone tissue/bio-material composite remodelling so facilitating the gradual substitution of the bio-material with natural tissue. Moreover, it is seen in which ranges of the aforementioned parameters residual regions of bio-materials remain in remodelled bone and when one can hope to get a final reconstructed bone constituted by pure bone tissue, i.e. when all bio-material is completely replaced by bone tissue. Many important problems are left open for further investigations in the theoretical development of the model, which seems promising for future researches. One of the most intriguing from the mathematical point of view concerns the description of bone tissue and bio-material micro-structure at different length scales by means of the introduction of higher gradient continuum models. However, the most urgent next investigations should concern i) the effects of the interplay of the different introduced parameters when their values are changed and the explanation of the different bio-mechanical phenomena described ii) the possibility of extending the proposed model by including the fluid flow phenomena occurring in the interconnected pores present in both bone tissue and bio-material. Indeed, evolution equations and natural boundary conditions which are able to describe fluid flow in porous media are already available in the literature (see e.g [START_REF] Dellisola | Boundary Conditions at Fluid-Permeable Interfaces in Porous Media: A Variational Approach[END_REF]) and can be suitably adapted to the present model. The introduction of a second or third dimension in the numerical examples may lead to new interesting observations as more complicated geometry would surely include more complex bio-mechanical synthesis and resorption mechanisms and scenarios.

  which is actually equivalent to equations (2) when considering a constant Poisson's ratio and when setting C b = C m = D b = D m = β. Indeed, it can be recognized that, in this case,

Figure 1 .

 1 Figure 1. Initial distribution of densities of bone and artificial material side (ρ * b /ρ max = 0.5 and ρ * m /ρ max = 0) and only artificial material on the right side ( (ρ * b/ρ max = 0 and ρ * m /ρ max = 0.5) ). We remark that, due to the fact that we are using a two-solid macroscopic mixture theory, at each point of the specimen (x axis) one needs to have both a material point of natural bone tissue and of artificial material, eventually with different densities. We start discussing the time evolution of mass densities distributions for low level of mechanical loading shown in Fig.2. In order to correctly interpret all the following figures, it must be remarked that the

Figure 2 .

 2 Figure 2. Time evolution of bone tissue and artificial material mass density for a low value of the applied external forces ( F ext = 1)

FinalFigure 3 .Figure 4 .

 34 Figure 3. Equilibrium configuration of bone tissue and artificial material mass density for a low value of the applied external forces ( F ext = 1)

FinalFigure 5 .

 5 Figure 5. Equilibrium configuration of bone tissue and artificial material mass density for higher values of the applied external forces ( F ext = 1.7)

Figure 6 .

 6 Figure 6. Equilibrium configurations for different values of the applied external forces

Table 1

 1 Values of the parameters used in numerical simulations

	ρ b /ρmax at t = 0	ρm/ρmax at t = 0	k	θ	sb	rb	rm
	0.5	0.5	4	1	0.1	0.1	0.15

Table 2

 2 Values of the parameters used in numerical simulations

	Bext	Ẽb	Ẽm	β	η	D	S * 0
	0	1	1	1.9	1	10	0.1