Next Article in Journal
Enhanced Antimould Action of Surface Modified Copper Oxide Nanoparticles with Phenylboronic Acid Surface Functionality
Next Article in Special Issue
Approaches for the Prediction of Leaf Wetness Duration with Machine Learning
Previous Article in Journal
Pipeline Inspection Tests Using a Biomimetic Robot
Previous Article in Special Issue
Discriminative Multi-Stream Postfilters Based on Deep Learning for Enhancing Statistical Parametric Speech Synthesis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Dipartimento di Scienze di Base ed Applicate per l’Ingegneria (SBAI), University of Rome La Sapienza, 00161 Roma, Italy
2
International Research Center for the Mathematics and Mechanics of Complex Systems (M&MoCS), University of L’Aquila, 67100 L’Aquila, Italy
3
Dipartimento di Ingegneria Meccanica e Aerospaziale (DIMA), University of Rome La Sapienza, 00184 Roma, Italy
4
Dipartimento di Ingegneria Civile, Edile-Architettura e Ambientale (DICEAA), University of L’Aquila, 67100 L’Aquila, Italy
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Submission received: 23 February 2021 / Revised: 6 March 2021 / Accepted: 8 March 2021 / Published: 10 March 2021
(This article belongs to the Special Issue Bioinspired Intelligence II)

Abstract

:
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.

1. 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]).

2. The Model Employed to Describe Bone Remodeling

2.1. 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 ,
  • 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
V v o i d = B * ϕ ( X ) d B *
Analogously, one can define the Eulerian porosity as
V v o i d = B n ( x ) d B
which is the corresponding current void density related to the current total space B . The two porosities, naturally are linked together by the relationship
ϕ ( X , t ) = n χ ( X , t ) J ( X , t )
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
E i j ( X , t ) = 1 2 u i , j + u j , i + u i , k u k , j
2.
the change of the Lagrangian porosity [32]
ζ ( X , t ) = ϕ ( X , t ) ϕ * ( X , t )
ϕ and ϕ * being the Lagrangian porosity corresponding to the current and the reference configuration, respectively. We can also express ζ by the equation
ζ = · [ φ ( U u ) ] = · w
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
ϕ * ( X , t ) = 1 ρ b * ( X , t ) ρ ^ b + ρ m * ( X , t ) ρ ^ m
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
E = E ˜ ( E , E , ζ , ζ )
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]:
E = 1 2 λ ( ρ b * , ρ m * ) E ii E jj + μ ( ρ b * , ρ m * ) E ij E ij + 4 α 1 ( ρ b * , ρ m * ) E ii , j E jk , k + α 2 ( ρ b * , ρ m * ) E ii , j E kk , j + 4 α 3 ( ρ b * , ρ m * ) E ij , i E kj , k + 2 α 4 ( ρ b * , ρ m * ) E ij , k E ij , k + 4 α 5 ( ρ b * , ρ m * ) E ij , k E ik , j + 1 2 K 1 ( ρ b * , ρ m * ) ζ 2 + 1 2 K 2 ζ , i ζ , i K 3 ( ρ b * , ρ m * ) ζ E ii
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 well-established 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
Y = Y b Max ρ b * ρ ^ b 2 + Y m Max ρ m * ρ ^ m 2
Y b Max and Y m Max being the maximal elastic moduli, the standard Lamé parameters can be expressed in terms of the apparent mass densities as
λ = ν Y ( ρ b * , ρ m * ) ( 1 + ν ) ( 1 2 ν ) , μ = Y ( ρ b * , ρ m * ) 2 ( 1 + ν ) ,
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
α 1 = α 2 = α 4 = Y ( ρ b * , ρ m * ) 2 , α 3 = 2 Y ( ρ b * , ρ m * ) 2 , α 5 = 1 2 Y ( ρ b * , ρ m * ) 2
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:
K 1 = 1 S p = ϕ * K f + ( α B ϕ * ) ( 1 α B ) K dr 1
where K f is the stiffness of the fluid inside the pores and
K dr = Y 3 ( 1 2 ν )
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
K 3 = α B K 1
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
α B = a 1 φ * + ( 1 a 1 )
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 deformation E ˙ as
2 D s = 2 μ v E ˙ ij E ˙ ij 1 3 E ˙ ii E ˙ jj + κ v E ˙ ii E ˙ jj
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
B * δ E d B * + B * D s E ˙ ij δ E ij d B * = B * δ W ext d B *
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:
δ W ext = τ B * τ i δ u i d S * + T B * T α δ u α , j n j d S * + Ξ B * Ξ δ ζ d S * ,
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
E int = 1 2 K u [ [ u ] ] · [ [ u ] ] + 1 2 K u [ [ ( u ) n ] ] · [ [ ( u ) n ] ] + 1 2 K ζ [ [ ζ ] ] 2
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.

2.2. 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]:
ρ b * t ( X , t ) = A b S H ϕ with 0 < ρ b * ρ ^ b ρ m * t ( X , t ) = A m S H ϕ with 0 < ρ m * ρ m * ( X , 0 )
The functions A b and A m are assumed to be piece-wise linear:
A b , m S = s b , m S for S 0 r b , m S for S < 0
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]:
P X , t = B a 1 E s ( X 0 , t ) + a 2 0 t D s ( X 0 , τ ) d τ ϖ ρ b * ( X 0 , t ) e X X 0 2 2 D 2 d X 0 B e X X 0 2 2 D 2 d X 0
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 weight-bearing 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
S ( X , t ) = P ( X , t ) P ref s for P ( X , t ) > P ref s 0 for P ref r P ( X , t ) P ref s P ( X , t ) P ref r for P ( X , t ) < P ref r
where the two thresholds P ref s and P ref r delimitate the width of a lazy zone, i.e., an interval where the osteoregulatory balance of the bone tissue is preserved.

3. 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 B G B 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
f 1 ( X 2 , t ) = X 2 w 1 2 F 0 + F 1 sin ( Ω t )
where w is the width of the specimen, the other amplitudes are set as: F 0 = 1.68 × 10 3 Y b Max , 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].

4. 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 B G B 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 B G B 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.
Figure 2 and Figure 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 bio-resorbable 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.
Figure 6 and Figure 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.
Figure 8 and Figure 9 exhibit the initial and final distributions of elastic energy density E . The system becomes stiffer since the stored energy decreases during the process. The external actions clearly induce bending of the sample with a localization of the energy close to the longest verges. The effect of the interface is less evident than that of the ζ field.
Figure 10 and Figure 11 plot the initial and final distributions of the dissipated energy density. On the contrary of the elastic energy density, the dissipated energy density increases as the simulation progresses. We also note a different distribution of it, particularly a significant localization all along the boundaries of the graft.
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.

5. 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.

Author Contributions

Conceptualization, D.S., A.M.B. and I.G.; methodology, D.S., A.M.B. and I.G.; software, D.S.; validation, D.S., A.M.B. and I.G.; formal analysis and investigation, D.S., A.M.B. and I.G.; writing—original draft preparation, D.S.; writing—review and editing, A.M.B. and I.G.; supervision, A.M.B. and I.G.; All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
RMSroot mean square
PDproportional derivative
B G B bone-graft-bone

References

  1. Cho, S.H.; Andersson, H.M.; White, S.R.; Sottos, N.R.; Braun, P.V. Polydimethylsiloxane-based self-healing materials. Adv. Mater. 2006, 18, 997–1000. [Google Scholar] [CrossRef]
  2. Toohey, K.S.; Sottos, N.R.; Lewis, J.A.; Moore, J.S.; White, S.R. Self-healing materials with microvascular networks. Nat. Mater. 2007, 6, 581–585. [Google Scholar] [CrossRef]
  3. Eremeyev, V.A.; Skrzat, A.; Vinakurava, A. Application of the micropolar theory to the strength analysis of bioceramic materials for bone reconstruction. Strength Mater. 2016, 48, 573–582. [Google Scholar] [CrossRef]
  4. Eremeyev, V.A.; Pietraszkiewicz, W. Material symmetry group and constitutive equations of micropolar anisotropic elastic solids. Math. Mech. Solids 2016, 21, 210–221. [Google Scholar] [CrossRef]
  5. Eremeyev, V.A.; Pietraszkiewicz, W. Material symmetry group of the non-linear polar-elastic continuum. Int. J. Solids Struct. 2012, 49, 1993–2005. [Google Scholar] [CrossRef] [Green Version]
  6. Madeo, A.; dell’Isola, F.; Darve, F. A continuum model for deformable, second gradient porous media partially saturated with compressible fluids. J. Mech. Phys. Solids 2013, 61, 2196–2211. [Google Scholar] [CrossRef] [Green Version]
  7. Rosi, G.; Placidi, L.; dell’Isola, F. “Fast” and “slow” pressure waves electrically induced by nonlinear coupling in Biot-type porous medium saturated by a nematic liquid crystal. Z. Angew. Math. Und Phys. 2017, 68, 51. [Google Scholar] [CrossRef] [Green Version]
  8. Alibert, J.J.; Seppecher, P.; dell’Isola, F. Truss modular beams with deformation energy depending on higher displacement gradients. Math. Mech. Solids 2003, 8, 51–73. [Google Scholar] [CrossRef]
  9. Pideri, C.; Seppecher, P. A second gradient material resulting from the homogenization of an heterogeneous linear elastic medium. Contin. Mech. Thermodyn. 1997, 9, 241–257. [Google Scholar] [CrossRef] [Green Version]
  10. Abdoul-Anziz, H.; Seppecher, P. Strain gradient and generalized continua obtained by homogenizing frame lattices. Math. Mech. Complex Syst. 2018, 6, 213–250. [Google Scholar] [CrossRef] [Green Version]
  11. George, D.; Allena, R.; Remond, Y. A multiphysics stimulus for continuum mechanics bone remodeling. Math. Mech. Complex Syst. 2018, 6, 307–319. [Google Scholar] [CrossRef]
  12. George, D.; Allena, R.; Bourzac, C.; Pallu, S.; Bensidhoum, M.; Portier, H.; Rémond, Y. A new comprehensive approach for bone remodeling under medium and high mechanical load based on cellular activity. Math. Mech. Complex Syst. 2020, 8, 287–306. [Google Scholar] [CrossRef]
  13. Hernandez-Rodriguez, Y.; Lekszycki, T. Novel description of bone remodelling including finite memory effect, stimulation and signalling mechanisms. Contin. Mech. Thermodyn. 2020, 1–13. [Google Scholar] [CrossRef] [Green Version]
  14. Giorgio, I.; dell’Isola, F.; Andreaus, U.; Alzahrani, F.; Hayat, T.; Lekszycki, T. On mechanically driven biological stimulus for bone remodeling as a diffusive phenomenon. Biomech. Model. Mechanobiol. 2019, 18, 1639–1663. [Google Scholar] [CrossRef] [PubMed]
  15. Gross, T.S.; Edwards, J.L.; Mcleod, K.J.; Rubin, C.T. Strain gradients correlate with sites of periosteal bone formation. J. Bone Miner. Res. 1997, 12, 982–988. [Google Scholar] [CrossRef] [PubMed]
  16. Judex, S.; Gross, T.S.; Zernicke, R.F. Strain gradients correlate with sites of exercise-induced bone-forming surfaces in the adult skeleton. J. Bone Miner. Res. 1997, 12, 1737–1745. [Google Scholar] [CrossRef] [PubMed]
  17. Barchiesi, E.; Spagnuolo, M.; Placidi, L. Mechanical metamaterials: A state of the art. Math. Mech. Solids 2019, 24, 212–234. [Google Scholar] [CrossRef]
  18. dell’Isola, F.; Seppecher, P.; Alibert, J.J.; Lekszycki, T.; Grygoruk, R.; Pawlikowski, M.; Steigmann, D.; Giorgio, I.; Andreaus, U.; Hild, F.; et al. Pantographic metamaterials: An example of mathematically driven design and of its technological challenges. Contin. Mech. Thermodyn. 2019, 31, 851–884. [Google Scholar] [CrossRef] [Green Version]
  19. Di Cosmo, F.; Laudato, M.; Spagnuolo, M. Acoustic metamaterials based on local resonances: Homogenization, optimization and applications. In Generalized Models and Non-Classical Approaches in Complex Materials 1; Springer: Berlin/Heidelberg, Germany, 2018; pp. 247–274. [Google Scholar]
  20. Spagnuolo, M. Circuit analogies in the search for new metamaterials: Phenomenology of a mechanical diode. In Nonlinear Wave Dynamics of Materials and Structures; Springer: Berlin/Heidelberg, Germany, 2020; pp. 411–422. [Google Scholar]
  21. Yildizdag, M.E.; Tran, C.A.; Barchiesi, E.; Spagnuolo, M.; dell’Isola, F.; Hild, F. A multi-disciplinary approach for mechanical metamaterial synthesis: A hierarchical modular multiscale cellular structure paradigm. In State of the Art and Future Trends in Material Modeling; Springer: Berlin/Heidelberg, Germany, 2019; pp. 485–505. [Google Scholar]
  22. dell’Isola, F.; Seppecher, P.; Spagnuolo, M.; Barchiesi, E.; Hild, F.; Lekszycki, T.; Giorgio, I.; Placidi, L.; Andreaus, U.; Hayat, T.; et al. Advances in pantographic structures: Design, manufacturing, models, experiments and image analyses. Contin. Mech. Thermodyn. 2019, 31, 1231–1282. [Google Scholar] [CrossRef] [Green Version]
  23. Vangelatos, Z.; Melissinaki, V.; Farsari, M.; Komvopoulos, K.; Grigoropoulos, C. Intertwined microlattices greatly enhance the performance of mechanical metamaterials. Math. Mech. Solids 2019, 24, 2636–2648. [Google Scholar] [CrossRef]
  24. Yildizdag, M.E.; Barchiesi, E.; dell’Isola, F. Three-point bending test of pantographic blocks: Numerical and experimental investigation. Math. Mech. Solids 2020, 25, 1965–1978. [Google Scholar] [CrossRef]
  25. Turco, E.; Misra, A.; Sarikaya, R.; Lekszycki, T. Quantitative analysis of deformation mechanisms in pantographic substructures: Experiments and modeling. Contin. Mech. Thermodyn. 2019, 31, 209–223. [Google Scholar] [CrossRef]
  26. Turco, E. How the properties of pantographic elementary lattices determine the properties of pantographic metamaterials. In New Achievements in Continuum Mechanics and Thermodynamics; Springer: Berlin/Heidelberg, Germany, 2019; pp. 489–506. [Google Scholar]
  27. Turco, E.; Misra, A.; Pawlikowski, M.; dell’Isola, F.; Hild, F. Enhanced Piola–Hencky discrete models for pantographic sheets with pivots without deformation energy: Numerics and experiments. Int. J. Solids Struct. 2018, 147, 94–109. [Google Scholar] [CrossRef] [Green Version]
  28. Eugster, S.; dell’Isola, F.; Steigmann, D. Continuum theory for mechanical metamaterials with a cubic lattice substructure. Math. Mech. Complex Syst. 2019, 7, 75–98. [Google Scholar] [CrossRef] [Green Version]
  29. Della Corte, A.; Giorgio, I.; Scerrato, D. A review of recent developments in mathematical modeling of bone remodeling. Proc. Inst. Mech. Eng. Part H J. Eng. Med. 2020, 234, 273–281. [Google Scholar] [CrossRef]
  30. Giorgio, I.; Spagnuolo, M.; Andreaus, U.; Scerrato, D.; Bersani, A.M. In-depth gaze at the astonishing mechanical behavior of bone: A review for designing bio-inspired hierarchical metamaterials. Math. Mech. Solids 2020. [Google Scholar] [CrossRef]
  31. Coussy, O. Poromechanics; John Wiley & Sons: Hoboken, NJ, USA, 2004. [Google Scholar]
  32. Biot, M.A. Mechanics of deformation and acoustic propagation in porous media. J. Appl. Phys. 1962, 33, 1482–1498. [Google Scholar] [CrossRef]
  33. Morgan, E.F.; Yeh, O.C.; Chang, W.C.; Keaveny, T.M. Nonlinear behavior of trabecular bone at small strains. J. Biomech. Eng. 2001, 123, 1–9. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Sansalone, V.; Martin, M.; Haïat, G.; Pivonka, P.; Lemaire, T. A new model of bone remodeling and turnover set up in the framework of generalized continuum mechanics. Math. Mech. Solids 2021. [Google Scholar] [CrossRef]
  35. Lekszycki, T.; dell’Isola, F. A mixture model with evolving mass densities for describing synthesis and resorption phenomena in bones reconstructed with bio-resorbable materials. Z. Angew. Math. Mech. 2012, 92, 426–444. [Google Scholar] [CrossRef] [Green Version]
  36. Lu, Y.; Lekszycki, T. Modeling of an initial stage of bone fracture healing. Contin. Mech. Thermodyn. 2015, 27, 851–859. [Google Scholar] [CrossRef] [Green Version]
  37. Lu, Y.; Lekszycki, T. Modelling of bone fracture healing: Influence of gap size and angiogenesis into bioresorbable bone substitute. Math. Mech. Solids 2017, 22, 1997–2010. [Google Scholar] [CrossRef]
  38. Bednarczyk, E.; Lekszycki, T. A novel mathematical model for growth of capillaries and nutrient supply with application to prediction of osteophyte onset. Z. Angew. Math. Und Phys. 2016, 67, 1–14. [Google Scholar] [CrossRef]
  39. Frasca, P.; Harper, R.; Katz, J.L. Strain and frequency dependence of shear storage modulus for human single osteons and cortical bone microsamples—Size and hydration effects. J. Biomech. 1981, 14, 679–690. [Google Scholar] [CrossRef]
  40. Madeo, A.; George, D.; Lekszycki, T.; Nierenberger, M.; Rémond, Y. A second gradient continuum model accounting for some effects of micro-structure on reconstructed bone remodelling. C R Mec. 2012, 340, 575–589. [Google Scholar] [CrossRef]
  41. dell’Isola, F.; Andreaus, U.; Placidi, L. At the origins and in the vanguard of peridynamics, non-local and higher-gradient continuum mechanics: An underestimated and still topical contribution of Gabrio Piola. Math. Mech. Solids 2015, 20, 887–928. [Google Scholar] [CrossRef]
  42. Auffray, N.; dell’Isola, F.; Eremeyev, V.A.; Madeo, A.; Rosi, G. Analytical continuum mechanics à la Hamilton–Piola least action principle for second gradient continua and capillary fluids. Math. Mech. Solids 2015, 20, 375–417. [Google Scholar] [CrossRef] [Green Version]
  43. Germain, P. The method of virtual power in the mechanics of continuous media, I: Second-gradient theory. Math. Mech. Complex Syst. 2020, 8, 153–190. [Google Scholar] [CrossRef]
  44. Epstein, M.; Smelser, R. An appreciation and discussion of Paul Germain’s “The method of virtual power in the mechanics of continuous media, I: Second-gradient theory”. Math. Mech. Complex Syst. 2020, 8, 191–199. [Google Scholar] [CrossRef]
  45. Cowin, S.C. Bone poroelasticity. J. Biomech. 1999, 32, 217–238. [Google Scholar] [CrossRef]
  46. Giorgio, I.; Andreaus, U.; dell’Isola, F.; Lekszycki, T. Viscous second gradient porous materials for bones reconstructed with bio-resorbable grafts. Extrem. Mech. Lett. 2017, 13, 141–147. [Google Scholar] [CrossRef]
  47. Mindlin, R.D. Micro-structure in linear elasticity. Arch. Ration. Mech. Anal. 1964, 16, 51–78. [Google Scholar] [CrossRef]
  48. Toupin, R.A. Elastic materials with couple-stresses. Arch. Rational Mech. Anal. 1962, 11, 385–414. [Google Scholar] [CrossRef] [Green Version]
  49. Allena, R.; Cluzel, C. Heterogeneous directions of orthotropy in three-dimensional structures: Finite element description based on diffusion equations. Math. Mech. Complex Syst. 2018, 6, 339–351. [Google Scholar] [CrossRef]
  50. Cluzel, C.; Allena, R. A general method for the determination of the local orthotropic directions of heterogeneous materials: Application to bone structures using μCT images. Math. Mech. Complex Syst. 2018, 6, 353–367. [Google Scholar] [CrossRef]
  51. Peng, L.; Bai, J.; Zeng, X.; Zhou, Y. Comparison of isotropic and orthotropic material property assignments on femoral finite element models under two loading conditions. Med. Eng. Phys. 2006, 28, 227–233. [Google Scholar] [CrossRef]
  52. Goda, I.; Assidi, M.; Belouettar, S.; Ganghoffer, J.F. A micropolar anisotropic constitutive model of cancellous bone from discrete homogenization. J. Mech. Behav. Biomed. Mater. 2012, 16, 87–108. [Google Scholar] [CrossRef]
  53. Alibert, J.J.; Della Corte, A. Homogenization of nonlinear inextensible pantographic structures by Γ-convergence. Math. Mech. Complex Syst. 2019, 7, 1–24. [Google Scholar] [CrossRef] [Green Version]
  54. Giorgio, I.; Harrison, P.; dell’Isola, F.; Alsayednoor, J.; Turco, E. Wrinkling in engineering fabrics: A comparison between two different comprehensive modelling approaches. Proc. R. Soc. Math. Phys. Eng. Sci. 2018, 474, 20180063. [Google Scholar] [CrossRef]
  55. Abali, B.E.; Wu, C.C.; Müller, W.H. An energy-based method to determine material constants in nonlinear rheology with applications. Contin. Mech. Thermodyn. 2016, 28, 1221–1246. [Google Scholar] [CrossRef]
  56. Rosi, G.; Placidi, L.; Auffray, N. On the validity range of strain-gradient elasticity: A mixed static-dynamic identification procedure. Eur. J. Mech. A/Solids 2018, 69, 179–191. [Google Scholar] [CrossRef] [Green Version]
  57. De Angelo, M.; Placidi, L.; Nejadsadeghi, N.; Misra, A. Non-standard Timoshenko beam model for chiral metamaterial: Identification of stiffness parameters. Mech. Res. Commun. 2020, 103, 103462. [Google Scholar] [CrossRef]
  58. Placidi, L. A variational approach for a nonlinear 1-dimensional second gradient continuum damage model. Continuum. Mech. Therm. 2015, 27, 623–638. [Google Scholar] [CrossRef]
  59. Placidi, L. A variational approach for a nonlinear one-dimensional damage-elasto-plastic second-gradient continuum model. Continuum. Mech. Therm. 2016, 28, 119–137. [Google Scholar] [CrossRef]
  60. Misra, A.; Singh, V. Micromechanical model for viscoelastic materials undergoing damage. Continuum. Mech. Therm. 2013, 25, 343–358. [Google Scholar] [CrossRef]
  61. Placidi, L.; Barchiesi, E.; Misra, A. A strain gradient variational approach to damage: A comparison with damage gradient models and numerical results. Math. Mech. Complex Syst. 2018, 6, 77–100. [Google Scholar] [CrossRef]
  62. Placidi, L.; Misra, A.; Barchiesi, E. Two-dimensional strain gradient damage modeling: A variational approach. Z. Angew. Math. Phys. 2018, 69, 1–19. [Google Scholar] [CrossRef]
  63. Placidi, L.; Barchiesi, E. Energy approach to brittle fracture in strain-gradient modelling. Proc. R. Soc. Math. Phys. A Eng. Sci. 2018, 474, 20170878. [Google Scholar] [CrossRef]
  64. Timofeev, D.; Barchiesi, E.; Misra, A.; Placidi, L. Hemivariational continuum approach for granular solids with damage-induced anisotropy evolution. Math. Mech. Solids 2020, 1081286520968149. [Google Scholar] [CrossRef]
  65. Giorgio, I.; Andreaus, U.; Scerrato, D.; dell’Isola, F. A visco-poroelastic model of functional adaptation in bones reconstructed with bio-resorbable materials. Biomech. Model. Mechanobiol. 2016, 15, 1325–1343. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  66. Garner, E.; Lakes, R.; Lee, T.; Swan, C.; Brand, R. Viscoelastic dissipation in compact bone: Implications for stress-induced fluid flow in bone. J. Biomech. Eng. 2000, 122, 166–172. [Google Scholar] [CrossRef] [Green Version]
  67. Cowin, S.C.; Nunziato, J.W. Linear elastic materials with voids. J. Elast. 1983, 13, 125–147. [Google Scholar] [CrossRef]
  68. Biot, M.A. Generalized theory of acoustic propagation in porous dissipative media. J. Acoust. Soc. Am. 1962, 34, 1254–1264. [Google Scholar] [CrossRef]
  69. Giorgio, I.; Andreaus, U.; Scerrato, D.; Braidotti, P. Modeling of a non-local stimulus for bone remodeling process under cyclic load: Application to a dental implant using a bioresorbable porous material. Math. Mech. Solids 2017, 22, 1790–1805. [Google Scholar] [CrossRef] [Green Version]
  70. Aretusi, G.; Ciallella, A. An Application of Coulomb-Friction Model to Predict Internal Dissipation in Concrete. In Mathematical Applications in Continuum and Structural Mechanics; Marmo, F., Sessa, S., Barchiesi, E., Spagnuolo, M., Eds.; Springer Nature Switzerland AG: Cham, Switzerland, 2021; Volume 127, Advanced Structured Materials. [Google Scholar]
  71. Heinemann, P.; Kasperski, M. Damping Induced by Walking and Running. Procedia Eng. 2017, 199, 2826–2831. [Google Scholar] [CrossRef]
  72. Eriksen, E.F. Cellular mechanisms of bone remodeling. Rev. Endocr. Metab. Disord. 2010, 11, 219–227. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  73. Green, A.E.; Rivlin, R.S. Multipolar continuum mechanics. Arch. Ration. Mech. Anal. 1964, 17, 113–147. [Google Scholar] [CrossRef]
  74. Polizzotto, C. A note on the higher order strain and stress tensors within deformation gradient elasticity theories: Physical interpretations and comparisons. Int. J. Solids Struct. 2016, 90, 116–121. [Google Scholar] [CrossRef]
  75. Andreaus, U.; Giorgio, I.; Madeo, A. Modeling of the interaction between bone tissue and resorbable biomaterial as linear elastic materials with voids. Z. Angew. Math. Phys. 2015, 66, 209–237. [Google Scholar] [CrossRef] [Green Version]
  76. Kumar, C.; Jasiuk, I.; Dantzig, J. Dissipation energy as a stimulus for cortical bone adaptation. J. Mech. Mater. Struct. 2011, 6, 303–319. [Google Scholar] [CrossRef] [Green Version]
  77. Andreaus, U.; Colloca, M.; Iacoviello, D. Optimal bone density distributions: Numerical analysis of the osteocyte spatial influence in bone remodeling. Comput. Methods Programs Biomed. 2014, 113, 80–91. [Google Scholar] [CrossRef]
  78. Andreaus, U.; Colloca, M.; Iacoviello, D. An optimal control procedure for bone adaptation under mechanical stimulus. Control. Eng. Pract. 2012, 20, 575–583. [Google Scholar] [CrossRef]
  79. Andreaus, U.; Colloca, M.; Iacoviello, D.; Pignataro, M. Optimal-tuning PID control of adaptive materials for structural efficiency. Struct. Multidiscip. Optim. 2011, 43, 43–59. [Google Scholar] [CrossRef]
  80. Carriero, A.; Pereira, A.; Wilson, A.J.; Castagno, S.; Javaheri, B.; Pitsillides, A.; Marenzana, M.; Shefelbine, S.J. Spatial relationship between bone formation and mechanical stimulus within cortical bone: Combining 3D fluorochrome mapping and poroelastic finite element modelling. Bone Rep. 2018, 8, 72–80. [Google Scholar] [CrossRef]
  81. Tiwari, A.K.; Prasad, J. Cortical Bone Adaptation to Mechanical Environment: Strain Energy Density Versus Fluid Motion. In Biomanufacturing; Springer: Berlin/Heidelberg, Germany, 2019; pp. 241–271. [Google Scholar]
  82. Hambli, R.; Kourta, A. A theory for internal bone remodeling based on interstitial fluid velocity stimulus function. Appl. Math. Model. 2015, 39, 3525–3534. [Google Scholar] [CrossRef]
  83. Forest, S. Micromorphic approach for gradient elasticity, viscoplasticity, and damage. J. Eng. Mech. 2009, 135, 117–131. [Google Scholar] [CrossRef]
  84. Cazzani, A.; Malagù, M.; Turco, E. Isogeometric analysis of plane-curved beams. Math. Mech. Solids 2016, 21, 562–577. [Google Scholar] [CrossRef]
  85. Cazzani, A.; Malagù, M.; Turco, E.; Stochino, F. Constitutive models for strongly curved beams in the frame of isogeometric analysis. Math. Mech. Solids 2016, 21, 183–209. [Google Scholar] [CrossRef]
  86. Greco, L.; Cuomo, M.; Contrafatto, L. A reconstructed local B formulation for isogeometric Kirchhoff–Love shells. Comput. Methods Appl. Mech. Eng. 2018, 332, 462–487. [Google Scholar] [CrossRef]
  87. Greco, L.; Cuomo, M.; Contrafatto, L. Two new triangular G1-conforming finite elements with cubic edge rotation for the analysis of Kirchhoff plates. Comput. Methods Appl. Mech. Eng. 2019, 356, 354–386. [Google Scholar] [CrossRef]
  88. Balobanov, V.; Niiranen, J. Locking-free variational formulations and isogeometric analysis for the Timoshenko beam models of strain gradient and classical elasticity. Comput. Methods Appl. Mech. Eng. 2018, 339, 137–159. [Google Scholar] [CrossRef]
  89. Yildizdag, M.E.; Ardic, I.T.; Demirtas, M.; Ergin, A. Hydroelastic vibration analysis of plates partially submerged in fluid with an isogeometric FE-BE approach. Ocean. Eng. 2019, 172, 316–329. [Google Scholar] [CrossRef]
  90. Yildizdag, M.E.; Demirtas, M.; Ergin, A. Multipatch discontinuous Galerkin isogeometric analysis of composite laminates. Contin. Mech. Thermodyn. 2020, 32, 607–620. [Google Scholar] [CrossRef]
  91. dell’Erba, R. Swarm robotics and complex behaviour of continuum material. Contin. Mech. Thermodyn. 2019, 31, 989–1014. [Google Scholar] [CrossRef] [Green Version]
  92. dell’Erba, R. Position-based dynamic of a particle system: A configurable algorithm to describe complex behaviour of continuum material starting from swarm robotics. Contin. Mech. Thermodyn. 2018, 30, 1069–1090. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Schematic representation of the system under study. It is constituted by two pieces of bone tissues joined by a graft made of a bio-resorbable material, that is the B G B junction.
Figure 1. Schematic representation of the system under study. It is constituted by two pieces of bone tissues joined by a graft made of a bio-resorbable material, that is the B G B junction.
Biomimetics 06 00018 g001
Figure 2. Plots of volume fraction of the bone tissue ρ b * / ρ ^ b at the end of the considered period for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Figure 2. Plots of volume fraction of the bone tissue ρ b * / ρ ^ b at the end of the considered period for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Biomimetics 06 00018 g002
Figure 3. Plots of volume fraction of the graft material ρ m * / ρ ^ m at the end of the considered period for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Figure 3. Plots of volume fraction of the graft material ρ m * / ρ ^ m at the end of the considered period for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Biomimetics 06 00018 g003
Figure 4. The initial distribution of the stimulus S as determined by the external actions for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Figure 4. The initial distribution of the stimulus S as determined by the external actions for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Biomimetics 06 00018 g004
Figure 5. The distribution of the stimulus S at the end of the considered period for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Figure 5. The distribution of the stimulus S at the end of the considered period for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Biomimetics 06 00018 g005
Figure 6. The distribution of ζ at the beginning of the time of simulation for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Figure 6. The distribution of ζ at the beginning of the time of simulation for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Biomimetics 06 00018 g006
Figure 7. The distribution of ζ at the end of the time of simulation for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Figure 7. The distribution of ζ at the end of the time of simulation for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Biomimetics 06 00018 g007
Figure 8. The distribution of the energy density E at the beginning of the time of simulation for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Figure 8. The distribution of the energy density E at the beginning of the time of simulation for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Biomimetics 06 00018 g008
Figure 9. The distribution of the energy density E at the end of the time of simulation for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Figure 9. The distribution of the energy density E at the end of the time of simulation for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Biomimetics 06 00018 g009
Figure 10. The distribution of the dissipated energy density at the beginning of the time of simulation for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Figure 10. The distribution of the dissipated energy density at the beginning of the time of simulation for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Biomimetics 06 00018 g010
Figure 11. The distribution of the dissipated energy density at the end of the time of simulation for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Figure 11. The distribution of the dissipated energy density at the end of the time of simulation for different substitute lengths of the graft: (a) 0.8 cm. (b) 0.6 cm. (c) 0.4 cm. (d) 0.2 cm.
Biomimetics 06 00018 g011
Table 1. Material parameters used in the numerical simulations.
Table 1. Material parameters used in the numerical simulations.
Y b Max (GPa) Y m Max (GPa) ν ρ ^ b (kg/m3) ρ ^ m (kg/m3)
1713.60.318001800
K f (GPa) K 2 (N) K u (N/m3) K u (N/m) K ζ (N/m)
1.7 1.7 × 10 5 1.7 × 10 18 1.7 × 10 7 1.7 × 10 7
κ v (N s/m2) μ v (N s/m2) s b (s/m2) r b (s/m2) r m (s/m2)
2.06 × 10 12 2.57 × 10 12 1.27 × 10 7 1.06 × 10 7 1.59 × 10 7
D (mm) a 1 a 2 P ref r (N/m2) P ref s (N/m2)
11150.9756.33
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Scerrato, D.; Bersani, A.M.; Giorgio, I. Bio-Inspired Design of a Porous Resorbable Scaffold for Bone Reconstruction: A Preliminary Study. Biomimetics 2021, 6, 18. https://0-doi-org.brum.beds.ac.uk/10.3390/biomimetics6010018

AMA Style

Scerrato D, Bersani AM, Giorgio I. Bio-Inspired Design of a Porous Resorbable Scaffold for Bone Reconstruction: A Preliminary Study. Biomimetics. 2021; 6(1):18. https://0-doi-org.brum.beds.ac.uk/10.3390/biomimetics6010018

Chicago/Turabian Style

Scerrato, Daria, Alberto Maria Bersani, and Ivan Giorgio. 2021. "Bio-Inspired Design of a Porous Resorbable Scaffold for Bone Reconstruction: A Preliminary Study" Biomimetics 6, no. 1: 18. https://0-doi-org.brum.beds.ac.uk/10.3390/biomimetics6010018

Article Metrics

Back to TopTop