Bio-Inspired Design of a Porous Resorbable Scaffold for Bone Reconstruction: A Preliminary Study

The study and imitation of the biological and mechanical systems present in nature and living beings always have been sources of inspiration for improving existent technologies and establishing new ones. Pursuing this line of thought, we consider an artificial graft typical in the bone reconstruction surgery with the same microstructure of the bone living tissue and examine the interaction between these two phases, namely bone and the graft material. Specifically, a visco-poroelastic second gradient model is adopted for the bone-graft composite system to describe it at a macroscopic level of observation. The second gradient formulation is employed to consider possibly size effects and as a macroscopic source of interstitial fluid flow, which is usually regarded as a key factor in bone remodeling. With the help of the proposed formulation and via a simple example, we show that the model can be used as a graft design tool. As a matter of fact, an optimization of the characteristics of the implant can be carried out by numerical investigations. In this paper, we observe that the size of the graft considerably influences the interaction between bone tissue and artificial bio-resorbable material and the possibility that the bone tissue might substitute more or less partially the foreign graft for better bone healing.


Introduction
Design of new materials with continuously increasing performance has always been a challenge for scientists and engineers. An approach rather promising is to study natural materials with non-common properties and try to reproduce their responses with new artificial solutions. In other words, the design of new synthetic materials is bio-inspired. For example, bone tissue, like other natural materials, is able to adapt its internal structure to the demands of the environment. Although some efforts to conceive a similar material have been made in recent past years, the way to cover for obtaining materials with the same features is still rather long [1,2]. The endeavor to conceive a material capable of mimicking the outstanding bone characteristic can be partially mitigated if one considers integrating the bone tissue with a synthetic graft. In this kind of 'composite' system, the primary objective is to design the graft in the best possible way to speed up the process of healing of an injured bone or, in the case of a disease that affects the tissue itself. On the one hand, some living bone tissue features will be inherited by the graft because of the system of cells that act already on that tissue and can, in some way, also perform on the graft material if this last is resorbable. On the other hand, it will be the inner microstructure of the graft that may accelerate or hinder the healing. Here, therefore, we focus on studying the interaction of the graft with the bone to provide guidelines for designing the structure of the implant. We already know that the graft must be porous and sufficiently strong to allow the bone remodeling and guarantee the mechanical strength required to carry on the entire process. However, in-depth knowledge of this kind of interaction is essential to improve the effectiveness of the graft to be designed. For all intents and purposes, that knowledge can be refined by employing suitable models of the system based on simplifying hypotheses that allow it to be described with sufficient accuracy and that will be corroborated experimentally. In the examined case, we have to deal not only with a complex mechanical system, such that the usual approach may be unsuccessful, but also with an evolution of the mechanical properties in the time due to the adaptation process taking place, and hence it is driven by mechanical and biological aspects. An efficient way to consider the typical complexity of the bone tissue, also with the addition of a graft intended to be as much as possible similar to the bone itself, is the use of generalized continuum theories (see, e.g., [3][4][5][6][7]) where the effects of the microstructure lying at the low levels of observation are introduced in the model as further kinematical descriptors that enhance the possibilities of the description of the system at the macro level. Together with poroelastic formulations, we believe that a second gradient theory should be considered because of the multi-scale organization of the bone tissue and its heterogeneity that involves the presence of high contrast in the mechanical properties of the constitutive elements of the microstructure [8][9][10]. This is a key factor to require the need of the second gradient continuum model. Besides, the consideration of the gradient of the deformation in the mechanical formulation is also beneficial for modeling the stimulus, which activates the adaptation process, [11][12][13][14] since it can be read as a macroscopic source of canalicular fluid flow that is commonly assumed to have a crucial role in the remodeling process [15,16].
To design the graft, a promising method is the one typically used to conceive metamaterials [17][18][19][20][21][22][23]. The opening gambit is the definition of the functional requirements of the graft, which are conveyed in a specific behavior to obtain. Then, the design effort is devoted to synthesizing the microstructure that can ensure the data specifications required. In this work, we assume that the graft has a microstructure very similar to the bone tissue; therefore, we model them, bone and graft, within the same formulation. However, we remark that the graft can be modeled with different mechanical properties to enhance its healing capability for future developments (see, e.g., [24][25][26][27][28]).

Mechanical Formulation
The system under study is constituted by a trabecular bone tissue to which an artificial graft, with a similar microstructure, is added to facilitate the healing process following a dire injury or a fracture. The main aim of the graft is to provide a junction that restores the mechanical continuity of the broken bone and a scaffold that allows the bone-cells responsible for the remodeling process to operate in resorbing the extraneous bio-material and substituting it with freshly produced bone. Since in this paper, we examine a preliminary case to illustrate as the proposed model can be beneficial from a computational point of view, bone is treated in general terms. Namely, no distinctions are made concerning the particular kind of bone (i.e., long bones, flat bones, etc.). The only hypothesis we consider is that the bone is trabecular because, in this case, the size-effect that a second gradient model can describe could be more evident. We are well aware that the effects of mechanical usage and biochemical agents probably occur at different rates and in various proportions for diverse skeleton regions. However, the hypotheses beneath this formulation are pervasive for all kinds of bone. From a modeling viewpoint, both the bone and the artificial graft can be described by a porous material. This is quite natural in the bone case, while some remarks should be made for the graft. Indeed, the bone tissue is characterized by different kinds of porosity, significant at various scales of observations, namely inter-trabecular, vascular, lacunar-canalicular, and collagen-apatite porosity [29,30]. In the case of the graft, instead, we can design it to the best following a suitable criterion that depends on the particular case examined. In most clinic cases, the objective is to use the graft as temporary support for the cells involved in the remodeling process that, at the end of the healing, will result in a more or less complete substitution of the graft with healthy bone. This perspective demands a graft as a material with inner cavities, recommendably interconnected, to supply the remodeling cells with the room to be able to carry out the process to which they are predisposed as well as with their nutrients.
Therefore, the investigation of the considered system is performed by using a generalized continuum model endowed by an extra-kinematical descriptor that takes into account the deformability of the micro-structure, namely the change of porosity, observed at the macro-level. In this framework, we assume that the significant behavior of the deformable body is completely determined when are known: the bulk displacement vector, u = x − X, 2.
the Lagrangian porosity, φ, where we used the standard definitions. Specifically, X represents any material particle of the body in the reference configuration (B * ), x = χ(X, t) is the place occupied by the particle X at the time t in the current configuration, and the Lagrangian porosity φ, assumed scalar, is the current void density referred to the reference total space, i.e., the current void space is evaluated as Analogously, one can define the Eulerian porosity as which is the corresponding current void density related to the current total space B. The two porosities, naturally are linked together by the relationship being J = det F = det ∇χ [31]. Therefore, with these kinematical descriptors the following deformation measures are allowed: 1.
the finite strain tensor, E, whose components are 2. the change of the Lagrangian porosity [32] ζ( φ and φ * being the Lagrangian porosity corresponding to the current and the reference configuration, respectively. We can also express ζ by the equation where U is the average fluid-displacement vector defined in such a way that the volume of fluid displaced through unit areas is ϕU and w = ϕ(U − u) represents the flow of the fluid relative to the solid but measured in terms of volume per unit area of the bulk medium.
Herein, we assume the finite strain tensor because of the microstructure that in some circumstances brings to a non-linear response nonetheless in small displacement and deformation regime as, e.g., it is documented in [33].
The considered system made of two parts, the bone tissue and the bio-resorbable graft, can be studied as a whole exploiting the mixture theory. Indeed, this theory has been introduced to describe a multiphase system made of several interpenetrable constituents. In other words, at any instant of time, all phases can be present at each material point. In particular, the bone tissue mainly has a solid phase and a porous space that is typically filled with bone marrow, interstitial fluids, blood, bone cells, etc. (see for more details [34]). The graft also is a porous material and when it is first placed onto the injured bone typically has empty pores. Those, soon enough, are filled with organic matter, like for the bone, in the healing process [35][36][37][38]. Therefore, after the first healing stage, the two parts are more or less in the same conditions. It is at this stage that we start to analyze the system. Clearly, even though the two considered parts start from a very similar initial configuration, their evolution due to the remodeling process is quite different. As a matter of fact, within the bone tissue, only a change of apparent mass density (the label apparent is used because the mass density is evaluated concerning the apparent or external volume of the body, including its pore volume) can be experienced. In contrast, in the bio-resorbable graft, together with the artificial material resorption (by osteoclasts), bone cells (osteoblasts) may synthesize new bone tissue. During the evolution of the remodeling process, inside the artificial graft, may coexist primarily three phases: bioresorbable material, bone, and physiological fluid insides pores.
In accordance with the mixture theory, the reference porosity can be expressed as follows depending on the volume fractions of the constituents, namely ρ * b /ρ b and ρ * m /ρ m , where ρ * b and ρ * m are the apparent mass densities of bone tissue and artificial material in the reference configuration, respectively; while the superimposed hat denotes the mass densities evaluated considering the true volume occupied by those phases. We remark that the apparent mass densities of the different phases evolve with time depending on the remodeling process and the mechanic environmental loads. Therefore, we assume the reference configuration of the continuum system is related to the apparent mass density in the reference state, i.e., when zero stress occurs. Eventually, the current porosity can be evaluated by A possible choice for the elastic energy density coherent with the assumptions made and frame-indifferent is The strain gradient tensor is adopted to model size and non-local effects observed in the considered continuum system [39][40][41][42][43][44]. The presence of ζ in the list of meaningful variables follows the suggestion by Biot [32] and Cowin [45] to model porous materials. The dependency on ∇ζ allows us to deal with boundary conditions related to the change of porosity. In classical poromechanics, this effect is usually neglected. Still, we believe that to model the interaction at the interface between bone and graft in detail, these kinds of boundary conditions should be taken into account. Indeed, a term in the energy involving ∇ζ permits the usage of the pore pressure at the bone/graft interface as a boundary condition between the two parts. In particular, we consider the following energy density [46][47][48]: where all the material parameters are considered a function of the apparent mass density of the bone and the graft evaluated in the reference configuration. Particularly, the coefficients α i are the second gradient rigidities and K j are the stiffnesses due to the porous nature of the system. Naturally, since they evolve with the remodeling process also the behavior of the considered system will change accordingly during the time. This is the simplest choice that one can conceive due to the isotropic assumption. Albeit, we are well aware that more complex models can be considered to describe the anisotropic nature of the bone [49,50], we refrain from considering this kind of complexity essentially because wellestablished remodeling process formulations are not yet available for these cases. Besides, the approximation due to the isotropic behavior, even though simplistic, in many real clinic instances, can be acceptable [51]. By introducing the Young modulus of the mixture as a power law of the solid-phases volume fractions and Y Max m being the maximal elastic moduli, the standard Lamé parameters can be expressed in terms of the apparent mass densities as assuming constant the Poisson ratio for the sake of simplicity. The material parameters related to the second gradient term of the energy density are estimated as follows Herein, the role of the characteristic length is of paramount importance to consider the effect of boundary layers due to the microstructure. The order of magnitude of is assumed to be 200 µm, i.e., the diameter of trabeculae that are the elements of which the microstructure is made of (see, for more details on lattice structures similar to the trabecular bone, [10]). The second gradient continuum model may be deduced using an homogenization technique that through a micro-macro identification permit to obtain an equivalent model that predicts at the macro scale the response of the material due to the microstructure existing at micro-scale (see e.g., [26,[52][53][54][55][56][57]). In this context, some critical issues related to the rise of damage in the system constituted by bone tissue and bio-material can be taken into account (see e.g., [58][59][60][61][62][63][64]).
The other three material coefficients K 1 , K 2 , and K 3 concern the presence of the pores and how they affect the macroscopic deformation equilibrium shapes of the system under study. Specifically, K 1 is a compressibility coefficient, namely the reciprocal of the poroelastic storage coefficient S p , which is defined as the volume of fluid released from unit bulk volume per unit decrease in pore pressure under the condition of constant confining stresses [32], and it can be written as: where K f is the stiffness of the fluid inside the pores and is the drained bulk modulus of the porous matrix, while α B (being φ * ≤ α B ≤ 1) is the Biot-Willis coefficient. K 2 represents the aptitude of the system to oppose to the establishment of a gradient of porosity and, for the sake of simplicity, is assumed constant. Finally, K 3 , as defined in [32], is expressed in the form being the coupling between the microstructure due to pores and the bulk solid. In the numerical simulations, we use a linear function of the porosity as with a 1 = 0.8. Albeit in the examined system, the sources of dissipation are many [65], foremost among those, we recall the dissipation: (1) in the bone solid matrix [66]; (2) in the fluid that fills the pores [67]; (3) at the interface between the solid and the fluid phase [68]; (4) at the interface between bone and graft [69]; (5) due to friction in possible micro-cracks within the solid phase [70]. Herein, we choose to treat all these sources from a macroscopic point of view, introducing a Rayleigh functional depending on the solid-matrix rate of deformatioṅ E as employing a Kelvin-Voigt model. In particular, in this formulation, the two contributions to the viscous dissipation, namely induced by a rate of shear deformation and a rate of hydrostatic deformation, are highlighted and particularized by the corresponding viscous coefficients µ v and κ v , respectively. It is worth noting that several time scales characterize the considered system. As a matter of fact, there are at least two characteristic times: the one specified by the application of external mechanical loads, which is about a few seconds (for example, in walking, the average frequency can be estimated below 2.0 Hz [71]); the cyclic period due to the remodeling process (for cortical bone it is about 120 days, while for trabecular bone it is about 200 days [72]). Naturally, the system is more complicated. One can also add the characteristic time due to the evolution of the stimulus that, once released, will remain in the bone tissue for some hours and then reabsorbed by metabolic processes. Therefore, a detailed formulation of the problem could be developed with a multiple time-scale technique, but the discussion of this approach is beyond the purview of this paper. Herein, we decide to consider only the slow time scale related to the remodeling process. The fast dynamic due to the mechanical loads is considered only for the effects induced at the slow time scale. To be more precise, the inertial effects are neglected, while a more thorough explanation of the other aspects is needed. Thus, to determine the mechanical load at the considered time scale, the root mean square (RMS) of the external force is considered. The RMS represents the value (say the amplitude) of a time-constant or slowly variable force that dissipates the same power of the original force, which is time-varying. The key idea is to consider the force at the slow time scale as it is commonly done in electrical engineering, defining the RMS current value equivalent to a direct current, which dissipates the same power in a resistor. In other words, we consider that the external forces have two time characteristic period, one fast and one slow: (1) the former is addressed attributing to the amplitude of the force its RMS value; (2) the latter is considered with a time dependence of the force at the slow time-scale. It is worth noting that the RMS evaluation is performed on a time interval equal to the period of the considered quantity if this is a sinusoidal or any periodic function. Otherwise, the time windowing of the RMS integration should be calibrated based on the significant frequency content of the signal. Namely, it should include a sufficiently high number of cycles of the smallest frequency of the interest. Similarly, in the Equation (17), the velocity amplitude should be considered as an RMS value, therefore the dissipation is taken into account for the fast scale resorting to the RMS operator while for the slow scale with the Rayleigh functional itself.
The mechanical model of the system is therefore obtained by applying the Generalized Principle of Virtual Work in the following form where δW ext is the virtual work done by external loads. Compatibly with the postulated energy density (10), we are allowed to consider these possible external actions: where the first term takes into account the forces per unit surface, i.e., τ i , on the part of the boundary ∂ τ B * , the second term is related to external double forces per unit of area [73], T α , which act on the normal part of the gradient of the displacement on the part of the boundary ∂ T B * [47,74], and the last term consider the effect of a pore pressure, Ξ, on the boundary ∂ Ξ B * . The external virtual works done by body forces per unit volume, like the weight, and forces per unit line on boundary edges are neglected.
To model the interactions between the bone tissue and the artificial bio-resorbable graft, an extra-contribution of the virtual work (18) is defined by integrating upon the interface surface the following potential where the symbol [[·]] represents the jump of a field f (X) from side to side through the interface, i.e., [[ f ]] = ( f + − f − ). This potential describes in particular an elastic junction between the two phases, bone and graft, that is established for the kinematical fields u, ∇u, and ζ. The material rigidities K u and K ∇u , and K ζ denote the effectiveness of the bond.

Growth/Resorption Process Formulation
This section details the model of the functional adaptation of the bone tissue with the inclusion of the bio-resorbable graft plugged in through surgery. Notably, we assume to describe this biological phenomenon by postulating the evolutive laws for the apparent mass density of the bone and graft in the reference configuration, namely when the stress within the system is zero. The simplest way to conceive these evolutions rules is by setting the time derivative of the mass densities as functions of a mechanical stimulus S(X, t) resulting from an external mechanical action and the current porosity φ(X, t). Typically, they are expressed as [35,75]: The functions A b and A m are assumed to be piece-wise linear: where the symbol {b,m} stands for the alternative between b, i.e., bone, and m, i.e., material of graft. The coefficients s and r are the rates of growth and resorption, respectively. Clearly, the artificial material of the graft cannot be synthesized by the bone cells (osteoblasts); hence the coefficient s m must be zero. The objective of the function H is to provide a tool for the calibration of limits of the apparent mass densities, namely a kind of end-of-stroke. In this context, the role of the porosity is vital. The pores allow the bone cells to reach the more remote zones of the system and operate to remove or synthesize bone tissue as demanded by the environmental external loads. Therefore, on the one hand, if the porosity is too low there is no room for this settlement, and the mass densities remain unchanged. On the other hand, if there is no physical support for the cells, i.e., high porosity, equally the remodeling cannot have a place. For these reasons, we can use a ∩-like shape function with its maximum approximatively near φ = 0.5 where the most effective conditions in the remodeling process can be supposed (see for more details [65]).
The stimulus S is the key variable to model the process of functional adaptation. It acts as a feedback correction to the mechanical properties of the system. In other words, the bone system is capable to monitoring its mechanical state and integrity. It is commonly accepted that osteocytes are supposed to perform this measuring. This information is used to compare the current condition with a reference state and a possible deviation between these two states triggers a physiological activity of the bone cells to restore the lost equilibrium. Herein, we choose that the osteoblasts and osteoclasts change the apparent mass densities of the solid phases and, in turns, the mechanical rigidities characteristics of the system.
The acquired information of the mechanical state by the osteocytes represents a detailed picture of the bone/graft system. From a modeling viewpoint, this knowledge can be expressed in the integral form [65,76]: which is a weighted average of the 'mechanical state' of the system. The exponential function defines the influence range of the osteocytes disseminated within the bone tissue, indeed. The radius of this area of influence is D. The farther the distance between the position of the osteocytes X 0 and the evaluation point of the signal X that will drive the osteoblasts and osteoclasts, the less powerful the signal level will be. The function denotes the density of osteocytes and acts as a gain for the 'acquired' mechanical state. Unquestionably, the more osteocytes there are, the better is for reconstruct the signal. Since the space density of the osteocytes is almost uniform in the bone tissue is reasonable to assume that is proportional to the bone density. In the numerical simulations, we set = 0.2 ρ * b /ρ b . In Equation (23), the mechanical state is computed as a linear combination, whose coefficients are denoted by a i , of the strain energy density E s of the solid matrix-including the first and second gradient displacement terms of Equation (10)-and the dissipated energy evaluated as the time integral of the dissipated power D s . This particular option is dictated by the main mechanical objectives of the bone system, namely to provide weightbearing capabilities, leverages for motion, and protection for the vital organs. Therefore, static strength can be assured looking to stored energy as a global figure of merit, while dynamic characteristics as shock-absorbing and vibration reduction could be achieved taking into account dissipative phenomena that can be linked to the rate of deformation. From a control point of view, this is very similar to a common proportional-derivative (PD) strategy (see, e.g., [77][78][79]).
We remark that the second gradient term in the present formulation is not only important for mechanical purposes but also to better describe the mechanism of functional adaptation. On the one hand, without this contribution we loose the possibilities of describes possibly boundary layers in the deformed shapes of the bone/graft system and size effects (the importance of these aspects is related to the characteristic length ); on the other hand, the strain gradient can be considered as a macroscopic source of canalicular fluid flow that seems to have a strong correlation to the formation of bone in specific sites (see, e.g., [15,16,[80][81][82]).
Finally, the stimulus can be written as (24) where the two thresholds P s ref and P r ref delimitate the width of a lazy zone, i.e., an interval where the osteoregulatory balance of the bone tissue is preserved.

An Illustrative Theoretical Case: Numerical Implementation
Although the proposed formulation is developed for the more general three-dimensional case, in the next sections, we consider a representative bi-dimensional example to better understand the potentialities of the proposed approach. Specifically, since the graft size typically varies from a few millimeters to some centimeters, we take into account a rectangular sample of size: 2 × 0.5 cm. We intend to study the interaction between the bone living tissue with a junction made with an artificial material that can be resorbed by the bone cells. This study aims to obtain beneficial criteria in the design of synthetic grafts to be used in bone reconstruction surgery. The sample comprises three sectors, two of bone tissue and a central one representing the graft, let say BGB junction for short. One side is constrained in order to avoid longitudinal displacement and any rigid motion keeping free the transverse displacement (see Figure 1). On the opposite side, a force per unit line with a symmetric linear distribution along the transverse direction is applied, as shown in Figure 1. The mechanical load is harmonically variable with a low frequency Ω, which directly affects the system at the considered time scale. We remind that the effects due to the fast oscillation of the force are incorporated in the model interpreting the magnitude of the force as an RMS value. In particular, we consider the longitudinal component of the force as follows where w is the width of the specimen, the other amplitudes are set as: F 0 = 1.68 × 10 −3 Y Max b , F 1 = F 0 /2 and Ω = 8.27 × 10 −6 Hz, which corresponds to five cycles per week. The transverse force component is assumed zero.
The material parameters used in the simulations are summarized in Table 1.   The numerical simulations have been performed through a finite element analysis using the COMSOL Multiphysics software which is able to solve the coupling problem of Equations (18) and (21). The discretization mesh size is set after a 'convergence' analysis to provide an adequate degree of accuracy. In particular, when the total energy of the sample remains almost unchanged, varying the size of the mesh, we accept that value, namely D/4. We also remark that a particular numerical technique has been used to properly consider the second gradient terms of the energy in the model. As commonly done for the beam formulation, a micromorphic approach based on the use of Lagrange multipliers allow us to treat the second gradient model in the same framework as the standard first gradient model [83]. Alternatively, one can implement an ad hoc numerical formulation that guarantees an inherent high continuity for the unknown fields, namely the isogeometric analysis (see, e.g., [84][85][86][87][88][89][90]). A novel numerical algorithm has recently been developed based on a particle system that has shown a beneficial capability to model second gradient materials and multi-crack evolution [91,92].

Results and Discussion
In order to illustrate one of the foremost parameters involved in the remodeling process and, therefore, the actual predictive capability of the proposed modeling tool to be used in the design of the graft, we carried out some numerical simulations. In what follows, we study the influence of the graft size inside the BGB junction, keeping the other parameters unchanged. In particular, we choose the length of graft: 0.8, 0.6, 0.4, and 0.2 cm. They correspond to the 40%, 30%, 20%, and 10% of the entire sample length.
The initial configuration is characterized by a uniform distribution of apparent mass density in the bone sides corresponding to a volume fraction of the bone tissue ρ * b /ρ b = 0.5; in the same way a uniform distribution of volume fraction for the graft material ρ * m /ρ m = 0.5 is fixed. Then, for the action of an external force, we promote bone growth in the BGB junction to analyze the effects of this mechanical stimulation by varying the size of the graft. The span of the time period analyzed is about 15 weeks to appreciate the evolution of the remodeling process. Figures 2 and 3 show the volume fractions of the bone and the graft material, respectively, at the end of the time interval of the simulations. The general trend we observe is that after an initial stage in which the artificial material is resorbed, at a different rate for the various cases in the graft part, the bone tissue starting from both lateral sides has also been synthesized in the central zone. Therefore, in the bone sectors, the volume fraction reaches saturation of the produced bone, and in the central zone, the two phases of bone and bio-resorbable material coexist at the same time. The different results in the behavior are a direct consequence of the initial distribution of the stimulus (see, Figure 4). Indeed, in the graft at the beginning of the process, there are not osteocytes, and hence the level of the stimulus may be negative depending on the size of the graft. In the cases where the length of the graft is smaller, the stimulus can be positive because it is a collective signal from all the osteocytes present in a certain space zone. Thus the osteocytes near the graft in the bone living tissue can produce a strong enough stimulus to be positive even in the graft zone. Naturally, if the graft is larger, the influence of the neighborhood osteocytes at the center of the graft fades. At the end of the considered period, instead, the distribution of the stimulus is almost the same for the four cases examined because the entire sample reaches a state of saturation, and the small differences that can be detected are related to the different amount of bioresorbable material which remains in the central zone (see, Figure 5). We notice herein that the artificial material and the bone tissue have a maximum (i.e., without pores) Young modulus slightly different. Besides, in the central zone where new bone tissue is also synthesized, new osteocytes are formed, and the stimulus level increases. At the end of the time interval of the simulation, notwithstanding the stimulus is positive, the remodeling process stops because the conditions for it are not any more favorable, i.e., the porosity is excessively small. Figures 6 and 7 display the initial and final distributions of the change of porosity ζ. From these pictures, it is clear that the process results in a stiffer system with a low level of porosity. Besides, we notice the differences due to the presence of the graft, especially at the bone-graft interface. The equilibrium shapes are plotted here and for the next figures, not in actual scale but with the same magnification (20×) to better recognize the kind of deformation. The numerical simulations show that for externally applied loads sufficiently high, the graft can be partially substituted by bone tissue. The composition of the final bone/graft composite in the graft zone primarily depends not only on the intensity of the mechanical load, which is a source of the stimulus, but also on the rates of bone synthesis and resorption of both phases and the size of the graft. As a matter of fact, the larger is the graft, the more resorption may occur far from the bone tissue where no osteocytes are present. This increase in porosity facilitates new bone production in the remodeling process that substitutes the graft material. Besides, the stimulus is nonlocal in nature [30], and therefore if the graft is sufficiently small, the influence due to the near living bone compensates for the absence of osteocytes within the graft at the beginning of the process. Thus, for small grafts, the synthesis of bone can start even before without or with a small extent of artificial material resorption. This qualitative behavior is in accordance with real bone healing cases, for which the process is favored for small gaps rather than large ones. In these cases, the nonlocal stimulus from the not damaged tissue triggers the healing faster and starts on the proximity of the healing process from which the stimulus originates. These results align with those of other works with similar formulations present in the literature (see, e.g., [35,40]).
We remark that the formulated model is based on the hypothesis that after the reconstruction surgery through the graft, this last is flood over by blood vessels, lymphatic vessels, nerves, and osteoblasts and osteoclasts responsible for the evolution of the system. After this stage, we have begun the simulations of the bone/graft evolution. To have a more detailed comprehension of the entire process, a proper description of the starting stage of the healing process is still missing in the proposed formulation. Therefore, this aspect is left open for further investigations.

Conclusions
This paper aims to show that with a relatively coarse model (without the details of the microstructure characterizing the system), it is possible to analyze the macroscopic bio-mechanical behavior of a system constituted by bone and artificial bio-resorbable material with decent accuracy. Therefore, after a brief introduction of the model employed, some numerical simulations are performed to confirm that the proposed formulation is qualitatively able to provide aid in the design of the graft to be used in bone reconstruction surgery. Naturally, the model parameters must be calibrated with experimental data, but this is always an obliged step. Assuming for the graft a mechanical micro-structure similar to that of the bone, we have shown that its size greatly influences the aptitude of the graft to be substituted with the newly synthesized bone. With this example in mind, we believe that the proposed model is a powerful tool for efficiently engineering the graft in an initial design phase by optimizing its macroscopic characteristics, which clearly requires many numerical simulations.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: RMS root mean square PD proportional derivative BGB bone-graft-bone