New description of gradual substitution of graft by bone tissue including biomechanical and structural effects, nutrients supply and consumption

A new description of graft substitution by bone tissue is proposed in this work. The studied domain is considered as a continuum model consisting of a mixture of the bone tissue and the graft material. Densities of both components evolve in time as a result of cellular activity and biodegradation. The proposed model focuses on the interaction between the bone cell activity, mechanical stimuli, nutrients supply and scaffold microstructure. Different combinations of degradation rate and stiffness of the graft material were examined by numerical simulation. It follows from the calculations that the degradation rate of the scaffold should be tuned to the synthesis/resorption rate of the tissue, which are dependent among the others on scaffold porosity changes. Simulation results imply potential criteria to choose proper bone substitute material in consideration of degradation rate, initial porosity and mechanical characteristics.


Introduction
Bone-substitute material has been intensively investigated over the past few decades. Traditionally, the bioinert materials such as stainless steel, titanium and its alloys have been used for fracture fixation and/or permanent grafts. However, there are some disadvantages of such materials. For fixation device, patient usually suffers from a secondary surgery for its removal after healing. For permanent graft, there exists a risk for unwanted bone resorption caused by the "stress shielding" due to the high stiffness of metallic material. Therefore, bioresorbable and/or biodegradable materials have recently attracted much attention for bone repair and tissue regeneration [6,36].
Graft bioresorbable material should satisfy several important requirements, among them: necessary strength and stiffness to work as mechanical support in early healing stage; promotion of tissue regeneration; gradual resorption and transfer of stress to newly formed bone tissue; and extra benefit of drug release accompanied with matrix resorption.
Among the most popular bioresorbable materials used for graft are polymers, ceramics and metals. Each of them has different inherent mechanical (strength and stiffness) and biological (degradation and cell-mediated resorption rate) properties. In addition, mechanical and biological performance of the scaffold is strongly affected by its microstructure described among the others by porosity, pore topology and inner surface area (total wall surface of the pores). The fact that microstructure determines porous material mechanical properties has long been understood, but the effect of biological influence is more complex. Microstructure of the scaffold is one of the most important factors affecting cell activities. On the other hand, cell activities contribute to the evolution of scaffold microstructure. The bone actor cells that synthesize (osteoblast) or destroy (osteoclast) the bone tissue are usually not present in the graft; they will need interconnected pore space to migrate into the scaffold and to proliferate there. Besides that, microstructure and porosity of surrounding material determine the ability of new blood vessels penetration into the pores. There exists a critical value of pore connections, above which the growth of blood vessels is possible [21,41]. The actor cells are functioning only on the matrix internal surface; therefore, the maximal possible actor cells number is determined by the inner surface area. Not only the cell number, but also the cell activity is directly controlled by geometric factors [5,31]. Inner surface curvature affects the shape of attached cells, which in turn regulates the synthesis and resorption rates.
In addition to cell number and cell activity, the phenomenon of degradation of the biomaterial is included in the formulation presented in this paper. Decomposition of the graft is the sum effect of cell-mediated resorption and self-degradation [34,35]. Improper design of the bioresorbable and/or biodegradable graft may lead to second fracture caused by too fast degradation, or remaining material inclusion caused by too slow degradation. In order to achieve the best healing result, scaffold's chemical composition, structural design and the type and amount of additives need to be carefully studied. It is not an easy task to test thousands of these combinations in experimental study. Therefore, a mathematical model which could predict the phenomena/scenarios happening in healing could be an alternative and cheaper solution.
In this work, we proposed a novel mathematical model that includes the basic interactions between mechanical, biological and structural properties of the bioresorbable scaffold and the living tissue. Numerical simulations were performed to examine the influence of the resorption rate on the healing process.

Kinematics description
The porous media considered in this study can be seen as a homogeneous material which is made of a mixture of bioresorbable porous bone tissue and porous bioresorbable and biodegradable material. Its kinetic behaviour depends on both its microstructural distribution [10] and proper definition of the regeneration kinetics [24,32] occurring within the structure. Bioresorption process is dependent on cell activities, while biodegradation process takes place without the aid of cells. In order to describe the kinematics of considered system, we introduce a reference configuration Ω. The points in Ω represent material particles which consist of mixture of both kinds of biomaterials.
-The displacement field u(x, t) is the first basic kinematic field introduced in this work. It describes the placement of considered body. This field is necessary to calculate dependent fields, as strain and stress, used to estimate mechanical stimulus controlling cell activities. -The next basic kinematic fields considered in our formulation are ρ b (x, t), ρ m (x, t). They describe the apparent densities of bone and the graft biomaterial, respectively. Their values change in time and space in response to a non-local mechanical stimulus.
-Concentration of nutrients c(x, t) is the last introduced basic field. It describes the amount of available nutrients in position x at instant t.

Mechanical property and state equation
The considered body is a porous mixture, where the solid phase is composed of the living bone tissue and the bone-substitute material. In the pores interstitial fluid, blood vessels and actor cells are present. The apparent densities of bone ρ b (x, t) and the bone-substitute material ρ m (x, t) are used to calculate the porosity of the mixture according to the following expression, where ρ b max and ρ m max are bulk densities of bone and biomaterial. Young's modulus of the mixture could be calculated by the power law model [20,22,26,27], where Y b andY m are bulk moduli of bone and biomaterial, respectively, and β b and β m are suitable constitutive coefficients.
For quasi-static analysis, the inertial term can be neglected. State equation and boundary conditions are represented as follows, In our model, porosity is taken into consideration indirectly by introduction of the apparent densities of the bone tissue and the graft material, Eq. 1. An alternative approach is possible, where instead of apparent densities, porosities of both materials are considered as independent variables (see [2,12,13,15]). In such situation, formula Eq. 6 introduced in the next section should be modified.

Cell population
Osteoblasts and osteoclasts are responsible for bone tissue synthesis and cell-mediated resorption, respectively. Both of the two kinds of cells can be active as remodellers only when attaching on the walls of the pores [33,37]. Though the ratio between them depends on different mechanical and biological conditions, their total number is only related with the available inner surface of the pores. This available surface depends on the porosity of the mixture and the topology of the pores. Specific surface is defined as the total surface of pores walls in unit volume of porous material.
Martin [29] and Fyhrie [9] investigated the relationship between porosity and specific surface area in cortical bone and woven bone and proposed, respectively, a fifth-order and a second-order dome-shaped polynomials to represent it. Taking into consideration the topology, especially the curvature of the pores surface [28,30,31,40], we proposed a Gaussian function instead of the dome-shaped one (see Fig. 1). This function reflects the space limitation due to small porosity, as well as the cell activation caused by the cell deformation as curvature increases. In general, the negative effect of small space available to cells at low porosity is partially reduced by the increased activity of cells attached to a large curvature substrate. However, this is still not clear to what extent the cell activities are dependent on substrate curvature. The advantage of this formulation is its simplicity and the fact that it approximately reflects the essential relations between porosity and surface of pores. However, more precise formulation is possible based on the approach often employed for foams description using Potts model (as, for example, in papers [16,25]). The normalized actor cell density is written as follows, where θ and ϕ 0 are parameters associated with the function shape. The left end of the Gaussian curve reflects the space limitation due to the nature porosity of cortical bone. In spite of the fact that even dense cortical bone has some degree of porosity, the pore space is not available for actor cells. This porosity is associated with the spaces occupied by Haversian and Volkmann's canals, lacuna cavities and canaliculi, around 5% generally. When osteoblasts become trapped in the bone matrix secreted by themselves, they transform into another phenotype osteocytes. The cavities where osteocytes stay are called lacunas. They are connected to each other by canaliculi. Haversian and Volkmann's canals are the longitudinal and transversal canals that contain blood vessels and nerves. Since these pores are occupied by either osteocytes or blood vessels and nerves, little space is available for the actor cells.
The right end of the Gaussian curve illustrates behaviours associated with the preferential tissue growth on concave surface compared with planar or convex ones [28,40]. At the pore size which is orders higher than the diameter of the actor cells, the surface is flat. When the curvature is comparable to cell size, cell expresses smaller anchor force [31] that in turn leads to weaker adhesiveness. Due to the traction force between cells, those cells with weaker adhesiveness are pulled up from the substrate. Therefore, extracellular matrix is formed beneath those cells.
The driving force of this phenomenon proposed by Rumpler [30] is explained by the minimization of the interfacial energy via the minimization of the interfacial area in thermodynamic perspective.
In contrast to actor cells, the sensor cells density is directly related with the apparent density of bone tissue. This is because the osteocytes in cortical bone are almost regularly distributed along the concentric cylindrical layers around the Haversian canals, and osteocytes in woven bone are generally evenly distributed in bone tissue. Although osteocyte densities in these two kinds of bones differ from each other by 2-4 times, in process of bone regeneration only the second case is considered [18]. Normalized sensor cells density is given by, where η is the parameter with value varying from 0 to 1 (0 ≤ η ≤ 1).

Nutrients transport and consumption
Nutrients concentration at selected place depends on environmental blood vessels density and distance from nearby vessels. Since in the regeneration phase there is a great need of nutrients from almost all regions, the assumption is made that the growth of blood vessels is chemical signal independent and has no preference of growth direction. Therefore, it is assumed that the growth of blood vessels simply follows the diffusion equation. In addition, the diffusion of nutrients from the capillaries to surrounding tissue is orders faster than the capillary growth rate; the time needed for the second process is negligible. In this context, equations from discrete particle models with stochastic dynamics are worth mentioning (see, for example, [4]). Thus, the total nutrients transport and consumption is described as follows, where ϕ defines porosity, D defines coefficient of diffusion, and R is consumption term. Percentages of nutrients consumption by actor cells and sensor cells present at x depend on their maximal consumption rates and cell numbers per unit volume, which are given by, where V a and V s are maximal nutrients consumption rates of actor cell and sensor cell, K a and K s are nutrients concentrations for actor cells and sensor cells at Nutrients consumptions by a single actor and sensor cells are described by the Michaelis-Menten kinetics [3], Total consumption by sensor cells and actor cells in unit volume is given by, Assuming that cell activity is proportional to the ratio of actual consumption to maximal consumption, thus activities of actor cell and sensor cell could be represented, respectively, by

Mechanical stimulus
Similar as in the previous works [22,23], the mechanical stimulus is measured by strain energy density SED. It is defined by the difference between the integral SED and the reference stimulus S 0 which represents the equilibrium state. The sensor cells respond to mechanical signal and release signalling molecule. Concentration of signalling molecules released by one sensor cell at one point decreases with the distance from the sensor cell. The mechanical signal at arbitrary point is then accounted by the integral of signalling molecules released by all adjacent sensor cells.
In the present work, mixture of the biomaterial and bone tissue is assumed to be linear elastic and isotropic material. The strain energy density is given by In addition, sensor cell density and activity are included in this formula. The stimulus is given by, where the exponential function e −a·r represents the signal decreasing with distance r from sensor cell. Parameter a is the reciprocal of influence distance. It is worth noting that a possible generalization of this stimulus can be done also adding a term related to the dissipation as proposed in [11,14,15]. Further extension of description of stimulus generation would be possible by consideration of local effects present in single cells, namely deformation of cell membrane. New achievements in description of lipid bilayers behaviours could be used to describe the cell deformation and get deeper insight into cell mechanobiology (see, for example, [38]). Inspired by the observations of bone remodelling process, according to which bone does not react to mechanical stimuli if the stimulus intensity is under certain threshold, we introduce the stimulus inactive zone (− l, l) in bone regeneration process (called often in the literature "lazy zone").

Mass density evolution
Finally, we propose the synthesis/resorption rate based on the following facts: osteoblast and osteoclast have different secretion rate; bone tissue can be absorbed or synthesized, while biomaterials undergo only resorption process; bone tissue and biomaterial respond differently to osteoclast activity (Fig. 2).
Synthesis/resorption rates for bone tissue are given by, where M b is amplification coefficient. Density evolution is: Besides cellular-mediated resorption, degradation in some biomaterials exists. Mechanisms of biodegradation are various from enzymatic degradation, hydrolysis and ion exchange, depending on type of material. The actual degradation rate varies in time, because it depends not only on material's inherent degradation rate but also on active factor concentration in surrounding fluid and exposed area (interface between material and surrounding interstitial fluid, that is the inner surface), etc. In this work, degradation is simplified as a constant rate process. It contributes to the materials resorption together with the cell-mediated resorption. Synthesis/Resorption rates for biomaterial are given by, where M m is amplification coefficient and A m0 is the degradation rate. Density evolution is:ρ

Numerical simulation
Relationships between different effects described in the previous sections are schematically illustrated in Fig. 3. Changes in the apparent density affect the strength of the matrix, which in turn affects the signal intensity sensed by sensor cell and finally regulates the actor cell activity. On the other hand, changes in apparent density cause the variation of matrix porosity, which control the actor cell attached on the inner surface. Total effect of these two regulating passways leads to the increase or decrease in cell-mediated synthesis/resorption rate. In order to verify the proposed mathematical model, a 2D model has been developed to simulate the creeping substitution of biomaterial by living bone tissue. Numerical simulation was performed using FEM software COMSOL Multiphysics v.5.2.®. The geometry model is similar to the fragment of rabit tibia investigated in experimental work performed by Iwasashi et al. [19]. A 4.5 mm × 7 mm rectangular defect was made in proximal tibia of rabbit and then fitted with a piece of rectangular prism of porous biomaterial. Unlike the       bio-inert material used in that work, the material in the proposed model is assumed to be bioresorbable and biodegradable. Figure 4a is a 2D representation of this model, which is the cut plane through the axial direction. The graft implanted in the left hand side column is the region of interest, and that in right hand side column is considered as the control group. Upper and lower rectangles represent very stiff elements with Young's modulus of 10 GPa, Poisson's ratio of 0.3. These two elements are introduced to simulate 3D system (shown in Fig. 4a) with Lower boundary of the model is totally constrained. A force F = 10 N is applied on the upper boundary and represents the effective mechanical loading of daily activity. The parameters used in simulations are listed in Tables 1, 2, 3, 4. Initially, the value of bone tissue density was assumed to be 1500 kg/m 3 . Different combinations of biomaterial mass densities and degradation rates were examined. The mass densities of the biomaterial were assumed to be 200, 800, 1400 and 2000 kg/m 3 . Thus, the porosities of bone and graft domains could be calculated by Eq. 1. The initial porosity for the bone is 0.25, and for the graft are 0.9, 0.6, 0.3 and 0.0.
In addition, degradation rates of 1 × 10 −3 , 1 × 10 −4 and 1 × 10 −5 kg/(m 3 s) were tested. The initial values of normalized nutrients concentration in the bone tissue and in the graft were assumed to be c = 1 mol/m 3 and c = 0 mol/m 3 . Nutrients diffuse into the graft area from the bone and the bone marrow cavity, as shown in Fig. 4a.
A quadrilateral mesh was generated over the whole domain, as shown in Fig. 4b. During the simulation, mass densities of both bone tissue and biomaterial change intensively in the domain of graft. Therefore, the mesh size similar to the influence distance 0.4 mm, that is the distance of about 15 osteocytes, was chosen. The bone mass density in domain of bone shaft will update in a moderate rate. For the sake of computing efficiency, the mesh size for the rest domains is set to be 0.7 mm.

Results and discussion
The parametric study was performed in order to determine the influence of the parameters associated with the graft. As mentioned in the previous section, three different degradation rates and four initial apparent mass densities were assigned to the graft. Each of the 12 studies was calculated for the process lasting 70 days. Total mass distributions at the end of the calculation are shown in Fig. 5. For the large degradation rate [1 × 10 −3 kg/(m 3 s)] group, only the model with weak material (200 kg/m 3 ) healed well (Fig. 5a). The remaining three groups (Fig. 5b-d) failed to heal, leaving a gap between the growing fronts. Size of the remaining gap grows as the initial mass density of the graft increases. That is because the weaker the material, the larger the mechanical stimulus, then the larger growth rate. The activity of the actor cells is reduced in case of strong graft and associated low stimulus. In such a case, the degradation rate of the graft is larger compared to the synthesis rate of the new bone. This leads to unsuccessful healing. Only if the growing front reaches the central area before the biomaterial is totally absorbed or degraded, the reconstructed bone could properly heal. Otherwise, there will be a gap remaining in the centre. In case of the presence of the weak section or even a gap in the centre of the reconstructed bone (left branch), the model will bend to left under compression. The bending causes a redistribution of strain energy density. The largest SED is no longer concentrated along the extension line of the bone shaft, but a little spreads into the bone marrow cavity. This effect is consistent with the observation in the three models with gap remaining, where the bending effect is more obvious. The newly formed bone for models with moderate (800 kg/m 3 ) (Fig. 5b), stronger (1400 kg/m 3 ) (Fig. 5c) and strongest graft (2000 kg/m 3 ) (Fig. 5d) grew towards the marrow cavity.
The combinations of moderate [1 × 10 −4 kg/(m 3 s)] or small degradation [1 × 10 −5 kg/(m 3 s)] rates with stronger or strongest graft also led to improper healing, as shown in Fig. 5h, k, l. The bridging parts between the two bone fragments do not reach the same mass density as the healthy bone tissue. Besides, the graft material is not fully absorbed due to large initial density and the small cell-mediated resorption rate. Figures 6 and 7 show the mass distribution of bone tissue and graft material, respectively. Even though bridging parts in the graft area in some models (Fig. 5g, j) have the same mass density as the healthy bone, they are still built of a mixture of bone tissue and graft material. As for models with large initial graft density (Fig. 7h, k, l), the mass of bridging parts is composed mostly of the remaining biomaterial.
In all of the models with weak graft material (200 kg/m 3 ) and the model with moderate graft and moderate degradation rate, a total substitution of the biomaterial by the living bone tissue at the end of the calculation can be observed.
If the external loading is in physiological range, two main simultaneous processes are responsible for bone healing, namely graft self-degradation and synthesis of a new tissue by actor cells activated by strong mechanical stimulus. To keep the balance between the rate of degradation and the rate of synthesis, a large enough self-degradation rate of the graft material is required. For the weak graft, the mechanical stimulus is high and results in big activity of the cells building new tissue. In such a case, the maximum allowed resorption rate is larger than that for strong graft. The situation is more complex when the graft material is strong with big density. The beginning of the process is associated with simultaneous self-degradation and resorption by the actor cells associated with low mechanical stimulus. When the density and the strength of the material decrease, the mechanical stimulus increases and activates the cells responsible for the synthesis of bone tissue. Since the mechanical stimulus is still relatively small, the required degradation rate should not excess the rate of synthesis.
However, an optimal healing procedure requires not only the good bridging quality at final healing stage under tested loading condition (static loading), but also: a sufficient strength during the whole healing process in order to exclude the potential risk of the secondary fracture in case of dynamic loading; short healing time during the whole healing procedure; and quick formation of the preliminary connection (osseoconnection) in the border area between the bone tissue and the graft.
In Fig. 8, the total mass density evolution at points in the centre of the graft for the models in 12 different conditions is presented. Total mass density in the models with large degradation rate (1 × 10 −3 kg/m 3 s) undergoes a short period of decrease. Next the total mass densities increase mildly or extensively corresponding to the different mechanical stimulus level and actor cell number, recovering to former level or reaching the same level as the healthy bone tissue.
The total mass densities in the models with weak graft grow quickly. The transition zones appear around 10, 11 and 32 days, respectively, for small, moderate and large degradation rate cases. The total mass densities in the models with moderate graft also grow very quickly. In addition, the transition zone for the moderate graft group appears earlier than the weak graft group. Figure 9 represents the bone mass density evolution of the points near the border between the bone and the graft. Group with moderate graft (800 kg/m 3 ) has the greatest initial growth rate because it provides a relative high porosity (0.6) and a relative low strength (2.5 GPa) that led to a sum effect of large growth rate. This group keeps its leading position during the first 8 days. The large portion of newly formed bone tissue on the border means a good connection quality. Figure 10 represents the biomaterial mass density evolution of the points in the centre of the graft. Biomaterial at this point in the four models with large degradation rate are fully absorbed at the end of the simulated procedure. Besides, moderate degradation rate enables a full resorption of biomaterial in the model with weak graft.
Both degradation and cell-mediated resorption contribute to the decrease in the mass density of the graft. Since the assumed degradation is a constant rate process, changing the curves slope is only associated with cellular behaviours. The curves representing four models with large degradation rate have practically bilinear shapes. Large decomposition rates, at the beginning of the healing process for these models, result from simultaneous mechanisms, namely cell-mediated resorption and self-degradation. After rapid changes, the process slows down, reflected by the second sections of the curves. This effect can be explained by weakening of the graft due to the rapid decomposition during the initial stages. It leads to the growth of mechanical signal. Increasing the mechanical signal results in the decrease in the osteoclasts activities, meaning that the contribution of cell-mediated resorption in total decomposition of material is small. In such situation, selfdegradation plays the main role in material changes (as shown in Fig. 2). The same effect could also be observed in all the other models, distinctly or indistinctly, with different lengths of section and slopes.
In the groups with stronger and the strongest grafts, a delayed cell-mediated resorption is observed. Three characteristic stages of density evolution can be observed in the plots for these models, Fig. 10. At the beginning, cell-mediated resorption is limited due to the small actor cell population. With the progress of the healing process, newly formed or enlarged cavities caused by material resorption provide more surface for actor cell attachment. This in turn enlarges the contribution of cell-mediated resorption reflected by larger slope of the curve. The activities of the cells and self-degradation lead to a weakening of the graft and a growth of mechanical stimulus. For the curves representing mass density evolution in the models of stronger graft with moderate and small degradation rates, as well as in the models of strongest graft with moderate degradation rate, a second significant slope change is observed. The resorption rates drop down to the smaller values.
Therefore, we conclude that, under these hypotheses, the best combination of initial graft density and degradation rate from the 12 studied models is moderate graft (800 kg/m 3 ) with moderate degradation rate [1 × 10 −4 kg/(m 3 s)]. This model enables a good interpretation of the occurring phenomena present during the healing process of bone reconstruction.

Conclusion
The theoretical model presented in the previous sections includes several important biomechanical and structural effects, such as the relation between porosity of bone and graft, number of active cells and their activities, mechanical stimulus, resorption and synthesis rates, nutrients transport and consumption. Most of these effects are linked to each other, complicating the understanding of the healing process. The effects of these relations can be observed in numerical simulations result.
The following new elements were included in the proposed model: -Changing the mechanical and structural properties of bone and bone-substitute material is described by the total effect of densities evolution of bone and bone-substitute material. -Biodegradation of the bone-substitute material is included in this model. -Nutrients consumption rates of sensor cell and actor cell are considered separately.
Unlike in the previous model [23], in which the switch between Michaelis-Menten kinetics and linear regime is introduced to describe the different consumption regimes in sufficient or insufficient nutrients supply cases, the consumption mechanism of the two cases is described here by of the admissible range kinematics. This formulation is more realistic because it enables description of the situation when even in case of insufficient amount of nutrients their residue remains. Consumption behaviours of both two kinds of cells follow the Michaelis-Menten kinetics. However, their preferable ranges of nutrients concentration and hypoxia tolerance are different. Therefore, their activities could be written in the same fashion, but with different parameters.
Based on this model, we tested different parameters associated with the graft material and concluded that the best combination is moderate initial graft density (800 kg/m 3 ) and moderate degradation rate [1 × 10 −4 kg/(m 3 s)].
However, the result is only an approximate description, accurate prediction requires verification of the parameters used in numerical simulation. One possible improvement of this model is to include the relationship between degradation rate of biomaterial and surface area, instead of a constant solution rate. With this new feature, the material inclusion could be predicted, even though the model did not include the microscale structure (single trabecula and pore).
Further developments may require an appropriate design of the microstructure of the graft to be transplanted in order to fulfil suitable healing features. Formulation based on the assumption about isotropic material can be generalized by the application of the approach used in the analysis of composite fibre materials. Since the graft can be seen as a scaffold on which the remodelling process develops, it may be very useful to optimize the structure of the graft taking into account the above-mentioned extensions so as to favour healing. Possible microstructures and second gradient effects to be used are, for example, those presented in [1,7,8,17,39].