An orthotropic continuum model with substructure evolution for describing bone remodeling: an interpretation of the primary mechanism behind Wolff’s law

We propose a variational approach that employs a generalized principle of virtual work to estimate both the mechanical response and the changes in living bone tissue during the remodeling process. This approach provides an explanation for the adaptive regulation of the bone substructure in the context of orthotropic material symmetry. We specifically focus upon the crucial gradual adjustment of bone tissue as a structural material that adapts its mechanical features, such as materials stiffnesses and microstructure, in response to the evolving loading conditions. We postulate that the evolution process relies on a feedback mechanism involving multiple stimulus signals. The mechanical and remodeling behavior of bone tissue is clearly a complex process that is difficult to describe within the framework of classical continuum theories. For this reason, a generalized continuum elastic theory is employed as a proper mathematical context for an adequate description of the examined phenomenon. To simplify the investigation, we considered a two-dimensional problem. Numerical simulations have been performed to illustrate bone evolution in a few significant cases: the bending of a rectangular cantilever plate and a three-point flexure test. The results are encouraging because they can replicate the optimization process observed in bone remodeling. The proposed model provides a likely distribution of stiffnesses and accurately represents the arrangement of trabeculae macroscopically described by the orthotropic symmetry directions, as supported by experimental evidence from the trajectorial theory.


Introduction
Wolff's law statement concerns the concept of "bone functional adaptation" to mechanical loadings.It is a well-known principle named after the German anatomist and surgeon Julius Wolff which states that the micromechano-morphological properties of bone are primarily, albeit not entirely, a consequence of bone response to mechanical deformation (Wolff 1892).More specifically, the law asserts that bone will remodel itself over time in response to the mechanical loads placed upon it, becoming denser and stronger in areas subjected to increased stress while weakened and more susceptible to fractures in areas subjected to less stress.Wolff's law has important implications for the prevention and treatment of conditions such as osteoporosis, as weight-bearing exercise and other forms of physical activity that can help maintain bone density and reduce the risk of fractures.Furthermore, specific dental treatments utilize Wolff's law to enhance the alignment of teeth within the dental arches (Cornelis et al. 2021).Bone remodeling in the vicinity of osteomyelitis (Masters et al. 2019) affected bone could also be influenced by the changed local loading environment both during disease development and after surgical treatment (Lamm et al. 2015).
This remodeling process is ensured by the synergic activity of three types of bone cells with different functions in the growth, maintenance, and repair of bone tissue: osteoclasts, osteoblasts, and osteocytes.Osteoclasts are large, multinucleated cells that dismantle mineral tissue through a process called resorption.They are responsible for removing old, damaged bone tissue and releasing calcium and other minerals into the bloodstream.Osteoblasts are cells that synthesize and deposit new bone tissue through a process called ossification.They secrete collagen and other proteins that form the framework for new bone tissue and also help to mineralize the bone matrix with hydroxyapatite, i.e., a calcium phosphate mineral.Osteoblasts are responsible for the growth and repair of bone tissue and play a key role in maintaining bone density and strength.Osteocytes are mature bone cells, deriving from the differentiation of osteoblasts, embedded within the mineralized matrix of bone tissue buried into tiny spaces called lacunae.They play a crucial role in regulating bone metabolism and responding to mechanical stresses placed on the bone that involves a process called mechanotransduction. Osteocytes are a kind of sensor cell connected to one another and to the bone surface by long, branching structures called dendritic processes housed in an intricate network of canals between the lacunae, called canaliculi, and filled with interstitial (periosteocytic) fluid.They communicate with the other bone cells to coordinate bone remodeling and repair.The activity of all these cells is essential for maintaining a balance between bone resorption and bone formation and for ensuring that bone tissue is constantly being renewed, replaced, and fit for its mechanical task.
The activity of bone cells results in the constitution and continuous shaping of two types of bone tissue that can be found in the human body: cortical and trabecular bone.They have different characteristics and serve different functions.Cortical bone, also known as compact bone, is the denser, more resistant outer layer of bone that makes up the shafts of long bones, such as the femur.It is composed of tightly packed layers of mineralized collagen fibers and is responsible for providing structural support and protection for the body and vital organs.They form cylindrical substructures called osteons.Cortical bone, due to its mechanical properties and space distribution, is relatively resistant to bending and twisting forces and is well-suited for weight-bearing activities.Trabecular bone, also known as spongy bone or cancellous bone, is a lighter, less dense type of bone tissue that is found at the ends of long bones and in the interior of short bones.It is composed of a network of bony struts called trabeculae that provide support for the bone and help distribute weight and stress.Trabecular bone is more compliant and can better withstand compression forces, making it an essential component of the body's shock-absorbing system, also considering that it is filled with bone marrow that has damping properties.Both types of bone tissue are important for maintaining overall bone health and function (Giorgio et al. 2021).
The bottom line is that the structure of bone tissue maintains a delicate balance between providing strength and support to the body while minimizing the cost of nutrients and required energy for the living cells operating in it (Bednarczyk andLekszycki 2016, 2022).By adjusting the distribution of bone mass and arranging its microstructure in an efficient composite structure, bone tissue is able to achieve this balance and maintain its function and integrity over time.Bone achieves this balance in response to the mechanical strain induced by external loads to which it is subjected.This means that areas subjected to more significant stress, such as the weight-bearing bones in the legs, will have more mass density than areas subjected to less stress, such as the bones in the fingers.This can be achieved, for example, through changes in the diameters and numbers of the trabeculae (Zhao et al. 2018).Besides, bone can also achieve an efficient composite structure by re-arranging its microstructure to maximize strength while minimizing weight, specifically by reorienting the pattern of the trabeculae depending on the direction of the external forces applied (Kivell 2016).
Usually, to explain the reorientation of the inner microarchitecture of bone tissue, the so-called trajectorial theory (Lanyon 1974) is employed, which assumes that the principal stresses are one of the main factors in determining the structure of trabecular bone tissue.More precisely, it is observed that the directions of the cancellous microstructure tend to be aligned to the eigenvectors of stress.Observed discrepancies in this alignment is mainly due to the different kinds of external loads that actually are applied to the bone tissue in different situations: each of these loads obviously determine different directions of principal stresses.Therefore, the effective trabecular orientation has to be considered as the weighted average result produced by the several mechanical loads applied during the functioning of the system "bone structure".
From a modeling viewpoint, bone remodeling can be regarded as a feedback phenomenon in the control theory framework (Frost 1987;Turner 1991).The stiffness tensor and material strength can be considered controlled quantities in this biological process.Thus, a control loop is established to regulate those variables at a setpoint, i.e., the so-called homeostatic state.To this end, the amount and the distribution of bone mass are changed to guarantee optimal functionality of the bone tissue.Indeed, based on the external request for mechanical resistance, the active bone cells, i.e., process actuators, can resorb or synthesize bone tissue based on the information provided by the osteocytes in charge of sensing the mechanical state of the tissue.In an effort to explain this biological process, an early attempt has been made introducing an evolutionary equation for apparent mass density (bone tissue is a porous material) that relied on a mechanical feedback signal (Beaupré et al. 1990;Mullender and Huiskes 1995).This is because bone mass density is directly linked to its stiffness; therefore, any changes in mass density would also alter the elastic modulus.Based on empirical evidence, this relationship between apparent mass density and stiffness has been introduced phenomenologically (Carter and Hayes 1976;Currey 1988).Thus, the change in the mass density also results in an adaptive evolution of the elastic modulus that optimizes the mechanical response of the bone.This approach is particularly suited for an isotropic behavior of the considered tissue since the apparent mass density is a scalar quantity and, hence, does not directly imply any mechanical properties linked with the reorientation of the trabecular microstructure.Therefore, the initial models incorporating mass density evolution have been primarily utilized for isotropic materials.
To take into account also the reorientation evolution, thus, it has been proposed a direct evolution of the stiffness tensor through a remodeling tensor that, in this context, is related to variables such as the apparent mass density and fabric tensor, which are intricately linked to the porosity and directionality of the trabeculae (Cowin et al. 1992;Doblaré and Garcıa 2002).These formulations can also be applied to orthotropic or anisotropic materials because of the tensorial nature of the stiffness tensor.However, some approaches are still based on mass density, even in orthotropic scenarios.Using phenomenological relationships, the distinct elastic moduli can be evaluated from the mass density.One potential limitation of this proposal is that the orthotropy directions are predetermined.In Sarikanat and Yildiz (2011), these directions are postulated to be aligned with the stress isolines following the trajectorial theory.However, it is worth noting that this perspective does not entail any form of evolution for the orthotropy directions.
The feedback signal is usually evaluated using a mechanical stimulus that conveys the information of the mechanical state considered pertinent for the particular adopted formulation.This mechanical state can be specified by strain energy density (Huiskes et al. 1987), effective stress (Carter et al. 1987;Beaupré et al. 1990), strain peak (Turner 1998), or damage variable (Prendergast and Taylor 1994;Hambli 2014), to name a few.There are two types of stimulus formulation: local and nonlocal.Local models apply the evolution law to each point of the body individually without directly considering the mechanical state of nearby points.Conversely, nonlocal models provide an evolution for each point in the body incorporating data from that point and other points in a specified finite surrounding volume of material.Nonlocal models are based on convolution integrals (Lekszycki and dell'Isola 2012;Kumar et al. 2011;George et al. 2018George et al. , 2019) ) or diffusion equations (Giorgio et al. 2019;Scerrato et al. 2022) and guarantee the possibility of describing the interaction between bone tissue and graft of bio-resorbable material or between healthy bone and necrotic tissue in which osteocytes could be dead.This kind of interaction is prevented from happening in local models because no stimulus can be generated where osteocytes are absent.On the contrary, with the nonlocal effect, the stimulus originating in surrounding areas can reach zones without osteocytes and trigger the remodeling process, as in fact occurs.
Bone tissue is characterized by a variegated architecture with multiple levels of complexity.Due to this, the conventional continuum elastic theory fails to accurately describe its intricate mechanics in different circumstances.For this reason, generalized continuum theories can be adopted to imitate the behavior of the real tissue more accurately and efficiently.The presence of high contrast in the mechanical properties of the bone may lead to the use of second gradient models that can represent a more rich behavior in terms of deformation modes and boundary conditions sustainable (Fedele 2022;dell'Isola et al. 2022;Solyaev et al. 2020;Vazic et al. 2021Vazic et al. , 2023;;Sarar et al. 2023;Carter et al. 1987;Andres et al. 2001).Strictly related models are those that deal with fiber-reinforced composites (Steigmann 2012;Franciosi et al. 2019;Spagnuolo 2022) since the mineralization of the tissue starts from a network of collagen fibers in the new woven bone.
Another aspect that is essential in a refined description of bone tissue is porosity.Models that take into account this aspect can explain phenomena involving inner pore pressure and possible dissipation due to the interaction between the solid matrix and viscous fluids, namely bone marrow, interstitial fluid, and blood, contained in it.The dissipation can be introduced using phenomenological laws, as can be found in (Jankowski et al. 2022;Sessa 2023), or more fundamentally, using Rayleigh functional (Giorgio et al. 2016).Due to the various types of porosity found in bone tissue, it is possible to utilize models with double or even higher porosity (see, e.g., De Cicco and De Angelis 2020;De Cicco 2022).Naturally, it is possible to use more complex models, such as poroelastic materials with fiber reinforcement and potential fluid inclusion, to depict the mechanical behavior of bones accurately (Tomic et al. 2014;Grillo et al. 2015Grillo et al. , 2018;;Cuomo et al. 2022;Gazzo et al. 2020).
The implementation of a variational formulation represents an elegant and effective approach to comprehending the relationship between the evolution of bone micro-architecture and the underlying physical laws that govern their behavior.This technique proves particularly advantageous in addressing complex problems, such as those encountered in studying bones and growth phenomena (Grillo and Di Stefano 2023;Grillo and Di Stefano 2023).This powerful approach can also be used to address challenging problems, such as those involving damage, that are intrinsically connected to the evolution of bone tissue (Placidi et al. 2018(Placidi et al. , 2018(Placidi et al. , 2019;;Timofeev et al. 2021;De Angelis 2007).Numerous authors interpret the remodeling of bone tissue as a complex optimization problem in which the local properties of mass and stiffness are regarded as conflicting objectives that must be balanced.From this perspective, once again, the variational tools provide a powerful arsenal at our disposal to obtain the optimized substructure underlying the bone tissue (see, as general reference Bendsoe and Sigmund 2003;Eschenauer and Olhoff 2001, and as a more specific application to the bone Nowak 2010Nowak , 2020)).
In this paper, we propose a possible way to incorporate the fundamental aspect of reorientating the trabecular microarchitecture of bone tissue in an elastic continuum scenario at a macroscopic level of observation.The material behavior is assumed to be orthotopic.The governing equation for the remodeling process is deduced from a variational point of view through the definition of a generalized virtual work that takes into account both the mechanical response and the evolution of the trabecular microstructure.The stimuli used as feedback signals to guide the process are defined on an energetic basis and have a nonlocal formulation that relays on diffusion equations.

Kinematics
As happens in standard continuum mechanics, the position x occupied by any particle X at the time t ∈ [t 0 , t f ] in the current configuration is given by the placement map , where E 3 is the Euclidean 3D affine space, where the material particles in the reference configuration are denoted by X in the domain B * .
Aiming to describe the bio-mechanical behavior of a bone tissue as a continuum deformable body having orthotropic material symmetry, in each material point and at every time instant, we introduce a material symmetry rotation Q(X, t).
Equivalently, one may consider the three time-variable orthogonal unit vectors where e i is an orthonormal basis in the space of dis- placements of E 3 .Moreover, we assume that the organism in which the bone tissue is incorporated has a bone tissue control system activated by suitable stimuli, which are generated as a response of the tissue mechanical state and which determine its production and adsorption.
Therefore, to model the bio-mechanical evolution of the bone tissue, we use the following list of kinematical descriptors: 1.The displacement field, denoted by u(X, t) = (X, t) − X, 2. The angles characterizing the rotation Q with respect to an Eulerian basis, 3. The set of material parameters, which characterize the orthotropic quadratic elastic energy, 4. The set of stimuli that are needed to guide the remodeling process.
The spirit here is to consider a generalized continuum theory where the overall evolution of the system is given by the standard descriptor, i.e., the placement or, equivalently, the displacement u , as well as further ones that specify both the mechanical properties of the tissue and the process of the bio-mechanical transduction guided on the informationbased stimuli production.
Let us denote by , the variables specifying the orientation A 1 (X, t) and A 2 (X, t) , by p(X, t) the entire set of elastic moduli, and, finally, by S = {S p i } the vector whose compo- nents represent the bio-mechanical stimuli, being each of them associated with each modulus p i .
In the present model, the deformation measure is given in terms of the deformation gradient tensor F = ∇ , and specifically by means of the linearized strain tensor Because of the complexity of the system, we presently confine our further development to the particular significant case in which a bi-dimensional domain B * is taken into account.The interesting results we believe to have obtained will motivate the development of a full 3D model.In the case of 2D evolutions, as only one angle is sufficient to characterize the orientation Q defining the evolutionary reori- entation of the trabecular substructure, we can proceed as follows.We introduce the macro-rotation related to the polar decomposition of F in which the rotation R is characterized by an angle in the bi-dimensional problem and U gives the deformation.Then, assuming that the microstructure associated with the trabeculae can be represented by the angle , which gives, e.g., the orientation of A 1 , we can introduce the relative angle that represents a measure of misalignment between the material directions of symmetry and the current orientation of the infinitesimal material cube at the same location due to the application of the external load.

Principle of virtual work
A generalized virtual-work principle is postulated for the independent virtual displacements u , its virtual gradient ∇u , the virtual change in the parameters , and p i as follows: ( (5) The first line of ( 5) is representative of the mechanical contribution given by the stored strain energy, namely w m , and the external virtual work done by the traction f per unit line and the density of double force since we consider second gradient materials (Madeo et al. 2012;Giorgio et al. 2017;Ganghoffer et al. 2019).The last two contributions are related to the mechano-biological evolution driven by the mechanical response.In them, we can interpret c  ̇ and as remodeling couples, where c and are constitutive parameters that relate to the time rate of the evolution of and the evolution of , respectively.In addition, c p i ṗi and −A(S p i ) could be thought as remodeling actions that expend a mechano-biological work on the material parameter p i .Similarly to the previous term, c p i is a constitutive param- eter linked to the time rate of the evolution of the elastic modulus p i and A is a biological potential associated with the evolution of p i .It is worth noting that the weak contribu- tion related to the remodeling actions in the last term of ( 5) produces an evolutionary term very similar to the one used in Eq. ( 5) in Giorgio et al. (2019).
The key idea of the model is to define a total bio-mechanical work as the sum of two contributions: 1) the first one able to characterize the mechanical response of bone and 2) the second one to determine the evolution of the inner architecture of the tissue considering both the orientation of the substructure and its stiffness.Since the complex phenomena involving the bone tissue develop on a multilevel time scale and the remodeling process is happening very slowly, herein, we neglect the rapid dynamic associated with the bone mechanics assuming a clear separation of time scale between the fastest variations related to the application of mechanical loads and the slow process of the remodeling.For this aim, we consider loads slowly variable in time, and thus, we neglect inertial phenomena.
To be more specific, we can introduce the energy (Cyron and Aydin 2017;DiCarlo and Quiligotti 2002): where w m is the standard orthotropic second gradient mechanical deformation energy and w r is the mechano-bio- logical term which captures biologically induced changes of the bone substructure and, thus, is representative of the stress-free configuration of the tissue over the remodeling process.
To describe the mechanical response of the bone tissue, we take as our starting point the constitutive equations for orthotropic materials.In this paper, we adopt, for the strain energy density w m , the Galilean invariant function: in which, for simplicity, we consider two additive terms: one related to a first gradient behavior and the other to a second (6) w(F, ∇F, , , p, S) = w m (F, ∇F, , p) + w r ( , p, S) gradient one.Particularly, the first gradient contribution can be specialized as where J 1 , J 2 , J 4 , J 5 , J 6 , and J 7 are the invariants character- izing the material symmetry in linear elasticity, that are defined as follows: Assuming small deformations, the dependence of the strain energy given by ( 8) could be assumed quadratic in the displacement field at leading order, and therefore, a function possessing material orthotropic symmetry is (Spencer 1982, 1984, Spencer and Soldatos 2007) in which K, , 1 , 2 , and 1 , 2 are material coefficients and 1 , 2 , and 3 are coupling coefficients.We remark that the first two terms in (10) are associated with the hydrostatic and deviatoric contribution to the isotropic part of deformation.
Here, all the material parameters may be inhomogeneous.
To complete the definition of the elastic stored energy, we assume also a second gradient contribution.The reason for which we postulate this further contribution to the stored energy must be attributed to the trabecular substructure of bone.Indeed, in the regions characterized by a very low mass density, we can assume that the trabeculae are very thin and long with weak connections, which may imply the need, at macro level, of second gradient models (Giorgio et al. 2017;Abdoul-Anziz and Seppecher 2018;Abdoul-Anziz et al. 2019;dell'Isola et al. 2022).This typically occurs in areas subject to minimal deformation (i.e., high porosity) or in the case of a disease like osteoporosis.In this special cases, a simple energy function which involves some second gradient effects with orthotropic symmetry is (Steigmann anddell'Isola 2015, Giorgio et al. 2016) where 1 and 2 are second gradient moduli and the notation "," indicates partial differentiation.Moreover, the second gradient contributions can be split into the tangential stretch gradient and the curvature parts considering the decomposition: (8) w 1 = ŵ1 (J 1 , J 2 , J 4 , J 5 , J 6 , J 7 ) (9) , i are the geodesic curvatures induced by the deformation on the trabecula (Shirani and Steigmann 2020) (see, e.g., Fig. 1).The vectors n i are orthogonal to a i .With the decom- position ( 12), the second gradient energy density is assumed to be being explicit the terms which refer to tangential stretch gradient effects and bending deformation.The same splitting procedure can be made for the stiffnesses, so we have s 1 and s 2 for the tangential stretch gradient stiffnesses and b 1 and b 2 for the bending stiffnesses.
The mechano-biologic energy density, w r , responsible for the remodeling is considered as the summation of several contributions.Particularly, one, w R , is related to the reori- entation of the orthotropic symmetry direction A 1 (which is associated with the directions of trabecular substructure), and all the others gathered in w s are associated with the evo- lution of the material stiffnesses characterizing the mechanical response.Thus, we have: By specifying the stored strain energy as (10) and (11), we can specialize p as {p 1 , p 2 , … , p 13 } = {K, , 1 , 2 , … , b 2 }.With the definition (4), we assume that where is a mechano-biological parameter associated with the level of energy that can be stored in such a process.( 12) Fig. 1 Material directions A i in the reference configuration aligned to the trabecula pattern and their images a i under the tensor F for the section of a femur It is worth mentioning that, in a more general context, the angle can be replaced with the rotation Q and the relative angle can be evaluated approximately by since the deformation tensor U tends to I in the small defor- mation regime. 1 Therefore, we have In the mechano-biological process, some energy dissipation occurs over time during the transition between different orientations; therefore, we define a Rayleigh dissipation function also to take into account this aspect, namely where c is a kind of "viscous" parameter.Therefore, a term c  ̇  can be introduced in the virtual-work formulation, which implies a non-conservative phenomenon.
To define the evolution of all the material stiffnesses and coupling coefficients, we use the contribution w s which depends on them.In particular, for each material parameter p i (with i = 1, … , 13 ), we can define the biologic potential energy as where A(S p i ) can be interpreted as a remodeling action which affects the material parameter p i and S p i is the stimu- lus associated to such an evolution (Giorgio et al. 2019).The remodeling action is assumed to be a linear function of the stimulus as follows: where p i represents a mechano-biological gain to be tuned with experimental evidence and S 0 p i is a reference stimulus introduced to set the homeostatic state.The difference between the actual stimulus S p i and the reference one gives feedback on how the system is far from homeostasis and, therefore, drives the evolution to restore this functional point of operation.
Finally, gathering all the contributions for the stiffnesses, we postulate Similarly to the reorientation evolution, a Rayleigh dissipation function is introduced also for the materials coefficients p i in the form: with c p i the related viscous coefficient.

Evolution equations for stimuli
In this section, we postulate, independently from the principle of virtual work, a set of evolution equations for biomechanical stimuli.
We consider as an open, and important, problem the postulation of a unique variational principle from which this last kind of evolution equations can be deduced as consequence (see, Barchiesi and Hamila 2022, for a first attempt, even if in a slightly different context).Here, we consider the stimulus field (as depending, in our simplified model, on the material particle and time) as a modeling tool in the description of the mechanically guided remodeling action occurring in bone tissues.
We postulate that the mechanical state, acquired by sensor cells, i.e., osteocytes, determines this stimulus.This stimulus is then the activator of the biological agents which remodel the bone tissue.As the principal path for the transmission of the information acquired by sensor cells, concerns the action of paracrine factors, which are chemical agents, released and spread via the diffusion, we postulate a diffusive model to describe the evolution of the stimulus field.
Therefore, for each material parameter p i , we postulate the evolution of the corresponding stimulus as follows: where d p i is a damping coefficient, p i the diffusion coeffi- cient assumed constant and thus isotropic for the sake of definiteness and tractability, r(U p i ) is a source term, and s(S p i ) is an absorption function.Following the work (Branecka et al. 2022), where the evolution of different material coefficients is determined by the energy contribution associated with the same coefficient, here we set the source term r as proportional to the portion of the energy density relative to the same coefficient.For example, to define this contribution related to the shear coefficient in (10), we set in the bi-dimensional formulation (21) 1 Indeed, we can write: and, thus, an approximation of the relative angle = − .and, for the sake of simplicity, the coefficient of proportionality is assumed to be 1.As usually done in the diffusion equation, the sink term s is defined proportional to the stimulus itself as follows: being R the absorption coefficient related to the metabolic activity associated with the stimulus.

Numerical results and discussion
The proposed model has been tested through numerical simulations in some illustrative cases.In particular, we consider a rectangular sample made of trabecular bone tissue.In the first test, we analyze a cantilever plate under a uniform distribution of a shear force density (t) on the side opposite the clamp (see Fig. 2).The second test, instead, is a three-point flexural one, where we employ the symmetry of the system to examine only half of the specimen (see Fig. 10).Since our formulation is variational, the numerical implementation of the tests has been conducted with the finite element method using the principle of virtual work (5) directly inside the commercial software COMSOL Multiphysics.We perform simulations in Comsol for testing the proposed approach because of its high flexibility and ease of use.Indeed, it is possible to implement the governing equation directly using the weak formulation.However, the general-purpose nature of the software does not allow for specific code optimization for the considered problem.To achieve this, one can consider implementing isogeometric analysis (Greco and Cuomo 2021;Greco et al. 2021;Torabi and Niiranen 2023) or an ad hoc discretization (Turco et al. 2016;Turco 2018Turco , 2022) ) in a subsequent phase of the development of the model.The simulations has been performed for a rectangular domain of 75 × 25 mm dimension, and the initial mate- rial stiffnesses are listed in Table 1.The second gradient stiffnesses, for the sake of simplicity, are set to be equal, namely The coefficients related to the evolutionary process for the reorientation of the material symmetry directions and the stiffnesses are reported in Tables 2 and 3. Herein, for the sake of simplicity and tractability, we postulate that the evolution of a few stiffnesses is more significant than the others, and then, we neglect the adaptive process of these last that changes their magnitude.Therefore, the only stiffnesses we consider to evolve are the shear modulus and the stretching  stiffnesses in the material symmetry directions, i.e., 1 and 2 .Finally, Table 4 shows the coefficients employed in the diffusion equations of the stimuli corresponding to the stiffnesses that are assumed to adapt to external mechanical conditions.The initial values of the stimuli are set to zero and Neumann boundary conditions specify that the stimuli flow vanishes to have insulation from the outside.In what follows, on the boundaries for which the boundary condition is not specified, the weak contribution or the dual is set to zero, as typically done for the adopted formulation.

A rectangular cantilever plate bending under a uniform load per unit length
In the first case, we perform a numerical simulation aiming at reproducing the behavior of a rectangular plate of bone tissue under the typical arrangement of a cantilever system.We set the displacement vector of the left short edge zero and apply a uniform shear force per unit length on the opposite edge variable in time.The initial transient of the force is designed to gradually reach a steady-state condition and avoid sudden variations of the displacement that where 0 is the magnitude of the force and T s the exten- sion of the transient.In particular, T s = 60, 480 s and 0 = 1.25 × 10 5 N/m .This value corresponds to a steady- state condition that is reached for the strain energy density w m in the considered case.
The results of this test show that the initially uniform distributed stiffnesses tend to increase their value where needed, i.e., where there is a localization of the contribution of energy associated with them, and decrease where the same energy density drops.Figures 3, 4, and 5 report this evolutionary change in the distribution of 1 , 2 , and , respectively, together with the current final deformation.Indeed, the localization of the apposition of new bone tissue and its resorption is placed in the expected regions.
The other significant aspect of the proposed model, which is the possibility of having an evolution of the material symmetry directions, is shown in three-time subsequent steps (see Fig. 6): 1) at the beginning, where they are assumed parallel to the sides of the rectangular specimen; 2) at an intermediate step; and at the steady-state condition.
To better understand the magnitude of this result, we consider the streamlines generated by these material symmetry directions at the end of the process and compare them with the isostatic lines related to the eigendirections of the strain tensor as shown in Fig. 7.
At first sight, they seem slightly different, but a more quantitative careful comparison demonstrates that they are ( 26) the same almost everywhere except for very tiny regions (the green ones), as shown in Fig. 8. Actually, Fig. 8 displays the difference between the angles related to the vectors A 1 and A 2 and those associated with the eigenvectors of the strain.This plot allows seeing that these two fields share the same angle or an angle of ∕2 in most of the domain.Therefore, since the two eigenvectors and the two fields of material symmetry must be orthogonal to each other, this implies that the isostatic and material symmetry directions coincide, and at least they exchange their order where the difference is ∕2 , accordingly with the trajectorial theory.The slight changes in the angles in the green regions of Fig. 8 are indeed related to the exchange between the two vector fields passing through areas with a different representation, namely a sort of transition zone.We can notice in Fig. 9 that passing through this region the same direction of orthotropy A i (in the picture A 1 is indicated with blue arrows and A 2 with red arrows) stops being parallel to one eigendi- rection of the strain tensor to align with the other, which is rotated of ∕2 relatively to the first, and simultaneously the same switch happens to the other direction of orthotropic symmetry.

A three-point flexure test
The second case examined is a three-point flexure test.The effects of a test fixture having pins interacting with the bone specimen during the test are simulated with an elastic barrier of potential that impedes the sample from overlapping with them.Specifically, we use a frictionless penalty method to mimic the contact between the pins and the bone tissue.The specimen is placed on two fixed supporting pins set apart distribution of 1 at steady state after the application of the load with a specific span length equal to 11∕6 L = 137.5 mm.A third loading pin in the middle of the test sample is used to deform the bone tissue gradually descending from above with a given time law.Here, we use the same expression as the previous test, Eq. ( 26), replacing the amplitude of the shear force with the assigned displacement to the loading pin u 0 = 2 mm.Due to the symmetry of the problem, we consider just half of the rectangular bone piece imposing symmetry conditions on the displacement at the middle section where the load is applied, namely u 1 = 0.
Figures 11, 12, and 13 exhibit the evolutionary change in the distribution of 1 , 2 , and , respectively, together with the current final deformation.As expected, analogously to the previous case, a distribution of the stiffnesses reflecting the actual distribution of the energy density of the considered contributions is achieved.
Figure 14 shows the evolutions of the material symmetry directions as the remodeling process progresses at three significant time instants.Their behavior is characterized by the alignment of them with the isostatic lines of the particular test considered, as can be clearly demonstrated by Fig. 15 where the streamlines generated by the unit vectors A 1 and A 2 , as well as the isostatic lines, are displayed.Once again, their difference is very small and limited to narrow areas with a switch in the order that yields a difference of ∕2 , thus, a matter of representa- tion, as it is quantified in Fig. 16.

Conclusions
Bone remodeling is the process by which bones continually renew themselves throughout a person's life.It involves the dissolution and absorption of old bone tissue by resorbing the mineralized matrix, breaking down the collagen fibers, and forming new bone tissue in a coordinated and controlled manner.Bone remodeling is essential for several reasons.It allows bones to repair themselves after injury or damage, helps to regulate the amount of calcium in the body, and plays a crucial role in maintaining bone strength and integrity.Besides, it also allows the bone to adapt to changes in stress and load over time since the new bone tissue is laid down in a highly organized manner and remodeled to meet the demands of the surrounding tissue.The process is influenced by a variety of factors, including hormones, diet, aging, and mainly physical activity.This process is really complex and we are far from completely understanding all its aspects.However, in this paper, we propose a new variational formulation that, by using the mechanical influence of the external load, is capable of explaining the adaptive process of rearranging the network of interconnected bone struts in the trabecular bone as well as their ability to carry on the load to which it is subject.The mechanical model is based on a generalized continuum theory and is characterized by an orthotopic material symmetry, which is entirely plausible for this material (Allena and Cluzel 2018;Cluzel and Allena 2018).The generalization involves some second gradient effects that could be relevant in specific contexts where the bone mass density is rarefied, and the lattice arrangement of the trabeculae is made of very thin and long struts with weak connections.Expressly, we assume a generalized principle of virtual work, assuming two key aspects, one related to the mechanical behavior of the bone tissue and the other concerning the evolutionary process that allows changing the orientations of the trabeculae as well as their thickness and length at the macroscopic level of observation of the entire sample.To illustrate the effectiveness of our proposal, two numerical tests have been performed.The results are promising and encouraging since they allow us to recover the optimization process typical of bone remodeling in the more general situation of orthotopic materials.Indeed, we are able to obtain not only a likely distribution of stiffnesses in the considered sample but also a reorientation of the material symmetry related to the actual arrangement of the trabeculae at a micro-level as postulated by the trajectorial theory for which we have experimental evidence.Tomek was more than a scholar; he was a mentor, a guide, and a source of inspiration for all of us.His dedication to imparting knowledge was unparalleled, and his enthusiasm for mechanics ignited a passion within each of us.Under his guidance, we delved into the complexities of bone mechanics, exploring the intricate workings of bone remodeling with his expert suggestions and invaluable advice.His commitment to excellence was evident in every discussion and every interaction we had with him.He possessed a remarkable ability to make even the most complex topics accessible and engaging.With his patient and compassionate approach, he encouraged us to push the boundaries of our understanding, nurturing our intellectual growth and fostering a deep appreciation for the field of mechanics.On behalf of all those whose lives he touched, we extend our heartfelt thanks to Tomasz Lekszycki.His influence will continue to resonate within us, shaping our careers and the world we live in.May his legacy inspire us to strive for greatness, and may his memory forever remain in our hearts.One of us (FdI) wants to further acknowledge a deep connection to Tomek with the following reminiscence: "I met Tomek in Warsaw at IPPT PAN in the winter of 2003: he was in the audience during my series of lectures on smart structures.Walking in the snow at minus twentyone Celsius, I discovered that he knew about epistemology, physics, mathematics, history, and variational principles applied to various problems in mechanics.The strength of Polish school of theoretical mechanics became clear to me when we started a long discussion about problems that I believed were very difficult and Tomek decided we could solve them together.He visited my department several times, obliged me to visit Warsaw often and participate in many biomechanics meetings.I, in turn, obliged him to come to Rome for an extended period secluding him in my house until he managed to write with me a complex code with which we could prove that our ideas about bone reconstruction were predictive of some interesting effects.I believe that we could prove that it is possible to design a bone scaffold which allows for bone reconstruction and complete resorption of bio-resorbable material.We became friends by calling each other and also asking for advice about the problems of life, which we tried to solve using Lagrangian methods.Tomek was revealed to be a skilled leader of a research group when he moved from IPPT PAN to Warsaw University of Technology, where he managed to start experimental activities guided by numerical simulations and predictive theory.He taught me how stubborn and proud and susceptible Polish person can be, but also that if a Polish is your friend, you can truly trust him.I regret that he left this world before we could fully prove our ideas were right.This paper is a small step in this direction.We wanted to formulate a variational principle for tissue growth, if, as he firmly believed, he is now somewhere looking at us, he will see that his ideas were not wasted.We are trying our best, Ciao Tomek!" Rest in peace, Tomek.

Fig. 3
Fig. 3 Bending test: distribution of 1 at steady state after the application of the load

Fig. 5 Fig. 6 Fig. 7 Fig. 8
Fig. 5 Bending test: distribution of at steady state after the application of the load

Fig. 9
Fig. 9 Bending test: Enlargement of the vectors A 1 (blue arrows) and A 2 (red arrows) in the transition zone

Fig. 12
Fig. 12 Three-point flexure test: distribution of 2 at steady state after the application of the load

Fig. 14
Fig. 14 Three-point flexure test: evolution of the material symmetry directions for the bending cantilever test: a at the beginning; b at an intermediate stage; c at the steady-state condition

Fig. 15
Fig. 15 Comparison between: a streamlines generated by the unit vectors A 1 and A 2 ; b isostatic lines for the three-point flexure test

Table 1
Initial material parameters (GPa) K