<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">\makeatletter
\@ifundefined{HCode}                  
{\documentclass[screen,CRMECA,Unicode,biblatex,published]{cedram}
\addbibresource{crmeca20240989.bib}
\newenvironment{noXML}{}{}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\sfrac#1#2{{#1}/{#2}}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\def\xtbody{\noalign{\relax}}
\def\tabnote#1{\vskip4pt\parbox{.51\linewidth}{#1}}
\def\xxtabnote#1{\vskip4pt\parbox{.65\linewidth}{#1}}
\def\mathbi#1{\boldsymbol{#1}}
\def\xmathop#1{\mathop{\text{#1}}}
\def\rmDelta{\Delta}
\usepackage{etoolbox}
\usepackage{mathrsfs}
\usepackage{upgreek}
\usepackage{graphicx}
\usepackage{amsfonts}
\usepackage{amssymb} 
\usepackage{bm} 
\usepackage{slashed}
\csdef{Seqnsplit}{\\}
\def\botline{\\\hline}
\newcounter{runlevel}
\usepackage{multirow} 
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\let\MakeYrStrItalic\relax
\def\jobid{crmeca20240989}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
\def\rmmu{\upmu}
\def\mn{\phantom{$-$}}
\def\refinput#1{}
\def\back#1{}
\newcommand{\bnabla}{\boldsymbol{\nabla}}
\newcommand{\bgamma}{\boldsymbol{\gamma}}
\newcommand{\bdelta}{\boldsymbol{\delta}}
\newenvironment{inftab}{}{}
\def\Custref#1#2{\hyperlink{#2}{#1}}
\def\og{\guillemotleft}
\def\fg{\guillemotright}
\makeatletter
\g@addto@macro{\UrlBreaks}{\UrlOrds}
\gappto{\UrlBreaks}{\UrlOrds}
\makeatother
\DOI{10.5802/crmeca.291}
\datereceived{2024-12-09}
\daterevised{2025-02-11}
\dateaccepted{2025-02-18}
\ItHasTeXPublished
}
{
\documentclass[crmeca]{article}
\usepackage[T1]{fontenc} 
\usepackage{upgreek}
\usepackage{slashed}
\def\CDRdoi{10.5802/crmeca.291}
\def\newline{\unskip\break}
\def\sfrac#1#2{{#1}/{#2}}
\def\xxtabnote#1{\tabnote{#1}}
\def\xmathop#1{\mathop{{#1}}}
\def\xtbody{\tbody}
\def\hypertarget#1#2{\label{#1}}
}

\usepackage{moreverb}

\dateposted{2025-03-10}
\begin{document}

\begin{noXML}

\CDRsetmeta{articletype}{research-article}

\title{Simulating bone healing with bio-resorbable scaffolds in a
three-dimensional system: insights into graft resorption and
integration}

\alttitle{Simulation de la cicatrisation osseuse avec des
\'echafaudages bio-r\'esorbables dans un syst\`eme tridimensionnel :
aper\c{c}us sur la r\'esorption et l'int\'egration des greffes}

\author{\firstname{Rachele} \lastname{Allena}\CDRorcid{0000-0002-8715-9174}\IsCorresp}
\address{Universit\'e C\^{o}te d'Azur, Laboratoire Jean Alexandre
Dieudonn\'e UMR CNRS 7351, Nice, France}
\address{Institut Universitaire de France, France}
\email[R. Allena]{rachele.allena@univ-cotedazur.fr}

\author{\firstname{Daria} \lastname{Scerrato}\CDRorcid{0000-0002-8746-0856}}
\address{Dipartimento di Ingegneria Meccanica e Aerospaziale (DIMA),
University of Rome La Sapienza, Italy}
\address{International Research Center for the Mathematics and
Mechanics of Complex Systems (M{\&amp;}MoCS), University of L'Aquila, Italy}

\author{\firstname{Alberto} \lastname{Bersani}\CDRorcid{0000-0001-8525-7363}}
\addressSameAs{3}{Dipartimento di Ingegneria Meccanica e Aerospaziale
(DIMA), University of Rome La Sapienza, Italy}
\addressSameAs{4}{International Research Center for the Mathematics and
Mechanics of Complex Systems (M{\&amp;}MoCS), University of L'Aquila, Italy}
\address{Gruppo Nazionale per la Fisica Matematica (GNFM) of the
Istituto Nazionale di Alta Matematica (INdAM, Italy), Italy}

\author{\firstname{Ivan}\nobreakauthor\lastname{Giorgio}\CDRorcid{0000-0002-0044-9188}}
\address{Dipartimento di Ingegneria Civile, Edile-Architettura e
Ambientale (DICEAA), University of L'Aquila, Italy}
\addressSameAs{3}{International Research Center for the Mathematics and
Mechanics of Complex Systems (M{\&amp;}MoCS), University of L'Aquila, Italy}
\addressSameAs{4}{Gruppo Nazionale per la Fisica Matematica (GNFM) of
the Istituto Nazionale di Alta Matematica (INdAM, Italy), Italy}

\begin{abstract}
This study introduces a three-dimensional (3D) model for investigating
the interactions between trabecular bone tissue and bio-resorbable
grafts, focusing on their role in bone remodelling processes. Using
principles of poroelasticity, the model captures the mechanical
behaviour of both bone and graft as porous continua, with their
interaction influenced by time-varying mechanical loads and diffusive
mechanical stimuli. The stimuli, originating from the strain energy
density, propagate through the system, triggering bone formation and
graft resorption across distant regions. Numerical simulations reveal
the critical impact of load frequency and intensity on remodelling
efficiency, with higher values promoting improved bone density and
graft integration. This 3D approach provides information on optimizing
scaffold design, offering valuable guidance to improve clinical
outcomes in bone repair and regeneration procedures.
\end{abstract}

\begin{altabstract}
Cette \'etude pr\'esente un mod\`ele tridimensionnel (3D) destin\'e
\`a explorer les interactions entre le tissu osseux trab\'eculaire et
les greffes bio-r\'esorbables, en se concentrant sur leur r\^{o}le dans
les processus de remodelage osseux. En s'appuyant sur les principes de
la poro\'elasticit\'e, le mod\`ele reproduit le comportement
m\'ecanique de l'os et de la greffe en tant que milieux poreux, dont
l'interaction est influenc\'ee par des charges m\'ecaniques variables
dans le temps et par des stimuli m\'ecaniques diffusifs. Ces stimuli,
issus de la densit\'e d'\'energie de d\'eformation, se propagent
\`a travers le syst\`eme, d\'eclenchant la formation osseuse et la
r\'esorption des greffes dans des r\'egions \'eloign\'ees. Les
simulations num\'eriques mettent en \'evidence l'impact d\'eterminant
de la fr\'equence et de l'intensit\'e des charges sur l'efficacit\'e du
remodelage, des valeurs plus \'elev\'ees favorisant une densit\'e
osseuse accrue et une meilleure int\'egration des greffes. Cette
approche 3D fournit des informations pr\'ecieuses pour optimiser la
conception des \'echafaudages, offrant ainsi des recommandations utiles
pour am\'eliorer les r\'esultats cliniques en mati\`ere de r\'eparation
et de r\'eg\'en\'eration osseuse.
\end{altabstract}

\keywords{\kwd{Bone remodelling}\kwd{Bone--graft
interaction}\kwd{Porous resorbable scaffold}\kwd{Fick's laws of
diffusion}\kwd{Bio-mechanical stimulus}}

\altkeywords{\kwd{Remodelage osseux}\kwd{Interaction
os--greffe}\kwd{\'Echafaudage poreux r\'esorbable}\kwd{Lois de diffusion
de Fick}\kwd{Stimulus biom\'ecanique}}

\shortrunauthors

\maketitle

\end{noXML}

\section{Introduction}\label{Introduction}

Bone remodelling is a lifelong process that maintains bone strength and
functionality by replacing old or damaged bone with new tissue,
ensuring structural integrity of the skeleton. This process is driven
by the orchestrated activities of osteoclasts, which resorb bone, and
osteoblasts, which form new bone. At the core of this process are
osteocytes, which act as mechanosensors, detecting mechanical stimuli
and coordinating the actions of osteoclasts and osteoblasts. The
stimuli sensed by osteocytes originate from the mechanical loads
applied to bones during daily activities, such as walking or
exercising. Bone tissue responds by adapting its architecture and
density to meet the mechanical demands placed upon it, a phenomenon
known as mechanotransduction~\cite{Frost2004,Turner}.

Mechanotransduction involves converting mechanical signals into
biochemical responses, triggering bone remodelling processes that
either strengthen the bone by increasing its density or weaken it by
resorbing less-needed areas. Dynamic loads, particularly those from
physical activity, are essential to stimulate bone remodelling. These
loads induce mechanical strains within the bone that are perceived by
osteocytes. In regions subjected to high mechanical strain, osteoblasts
are activated to synthesize new bone, leading to localized bone
formation. Conversely, osteoclasts are activated to resorb bone in
regions of low mechanical strain, resulting in bone loss
\cite{Frost2004}. Thus, bone architecture and mass distribution are
optimized to withstand mechanical forces, preventing fracture or
deformation~\cite{Turner,Robling}.

The process of bone remodelling is influenced not only by mechanical
forces but also by biochemical and cellular factors. Calcium
homeostasis and repair of micro-damages from daily stresses are other
critical functions of bone remodelling~\cite{Frost2004,Lanyon}.
Furthermore, the complexity of this process is reflected in the
interactions between cells and their surrounding environment (see,
e.g.,~\cite{morozova2020mmcs} for a similar mechanism of functioning).
Various signalling pathways are involved, particularly the diffusion of
mechanical stimuli within the bone matrix. These stimuli guide
osteocytes in assessing the mechanical environment and regulating bone
adaptation accordingly~\cite{Cowin1991,Allena2023}.

The role of mechanical stimuli in bone remodelling has been the focus
of numerous theoretical and experimental studies. Traditional models of
bone remodelling often treated the mechanical stimulus as localized,
affecting only the region where the strain occurred
\cite{Mullender,Huiskes}. However, recent advances suggest that a
diffusive model better reflects the behaviour of actual bone tissue
\cite{giorgio2019bmm}. By considering the stimulus as a diffusive
agent, some models account for the influence of mechanical loads on a
larger area, improving the accuracy of predictions regarding bone
adaptation. The stimulus is assumed to propagate through the bone
matrix, influenced by the microstructure of the tissue, such as its
porosity. This propagation affects osteocyte activity and,
consequently, bone resorption and formation rates~\cite{cowin1999bm}.
These models consider the stimulus as a function of the strain energy
density generated by mechanical loads. In this context, the strain
energy density serves as the source of the stimulus, which diffuses
through the bone tissue, initiating the remodelling response in distant
regions. The diffusive nature of the stimulus allows it to propagate
throughout the tissue, allowing a coordinated cellular response
throughout the bone structure, even in areas that do not experience
direct mechanical loads~\cite{Allena2023,Mullender}.

In addition to biological factors, the structural characteristics of
bone tissue, such as porosity, play a significant role in bone
remodelling. Bone is a composite material composed of a solid phase
consisting of collagen and mineral components, and a porous phase
filled with fluids such as blood and bone marrow. The porous nature of
bone allows the diffusion of nutrients and signalling molecules, which
are essential for cell migration and function. Porosity also affects
the mechanical properties of bone, influencing how mechanical loads are
distributed within tissue~\cite{Allena2024,Pivonka}.

The remodelling process is crucial not only for maintaining bone health
but also for healing bone injuries and defects. In clinical settings,
the use of bio-resorbable scaffolds to repair bone defects has become
an area of significant interest. These scaffolds are designed to
integrate with bone tissue and support mechanical loads during healing.
For effective bone healing, the scaffold must have suitable mechanical
properties, such as stiffness and toughness, and must also be porous
enough to allow cell migration and diffusion of signalling molecules
\cite{Allena2023,Carter}. The porosity of the scaffold is particularly
important because it facilitates the invasion of osteogenic cells and
the formation of new bone while at the same time providing structural
support until the scaffold is resorbed and replaced by natural bone
tissue~\cite{lekszycki2012zm}. The porosity concurs to define the
microstructure of bone tissue and can also be replicated in the graft
to facilitate osteointegration~\cite{sheidaei2019mms}, but to have a
complete picture of the mechanical behaviour of a system treated with a
bio-resorbable scaffold, it could be necessary, under some conditions,
to generalize the mechanical description in the sense of
\cite{germain2020mmcs}, with generalized models that include micropolar
effects, grain interactions, inclusions, fibres, etc.\ 
\cite{chroscielewski2020ijss,eremeyev2020zamp,placidi2022ijss,dellIsola2015mms,federico2019mmcs,tepedino2023mmcs,casalotti2024zamp,barchiesi2024mmcs}. 
Biological systems, such as our system, composed of bone tissue, an
implant, cellular networks, and equipped by physiological processes,
exhibit highly complex and interconnected dynamics. These systems are
often influenced by a multitude of factors, including genetic,
environmental, and stochastic variations, which introduce inherent
uncertainties. Accurate modelling of such systems requires accounting
for these uncertainties, as they can significantly impact system
behaviour and predictions. By incorporating probabilistic approaches
and sensitivity analyses, it is possible to capture the variability and
unpredictability that are fundamental to understanding and effectively
simulating biological complexity (see,
e.g.,~\cite{soize2023mmcs,lavalle2023zamm,lavalle2024zamp} regarding
the role of uncertainties in multiscale mechanics).

Mathematical models have been developed to study the interaction
between bone tissue and bio-resorbable scaffolds. These models aim to
optimize the design of scaffolds by predicting how their mechanical and
geometric properties influence the bone healing
process~\cite{dissaux2022cppj,eremeyev2016sm}. For instance, a recent
study employed a diffusive model to analyze the interactions between
bone and a resorbable graft, focusing on how mechanical stimuli
influence the rates of bone formation and graft
resorption~\cite{scerrato2022s,stoltz2018jci}. The results of this
study highlighted the importance of scaffold porosity and mechanical
properties in promoting effective bone regeneration. By simulating
different scaffold designs, such models provide valuable insight that
can guide the development of more efficient bone repair strategies.

The present study represents a significant extension of our previous
work by introducing a three-dimensional (3D) system that includes both
bone and bio-resorbable grafts, enhancing the realism of the model in
simulating the bone remodelling process. The new model simulates the
interaction between bone tissue and the graft, considering how
mechanical loads affect both components over time. As in
\cite{Allena2024}, we use a diffusive signal to represent the
mechanical stimulus, which originates from the strain energy density
generated by dynamic loads and propagates through both the bone and the
graft. {Additionally, our model is developed within a
poroelastic framework. As such, it implicitly takes into account the
interplay between interstitial fluid flow and solid deformation. This
constitutes an important feature, as osteocytes are known to be
sensitive to fluid-induced shear stress. The pressure gradients that
result inside the porous bone matrix contribute to fluid flow, which
can cause shear stress on osteocytes and their dendritic processes.
Although our model does not explicitly resolve the flow of the
microscale canalicular network, it captures a global description of
bone remodelling, combining indirect fluid-driven stimuli and direct
solid-phase strains into a single mathematical framework. We account
for these effects through the poroelastic formulation, which inherently
couples the solid and fluid phases of bone and graft materials. The
poroelastic framework considers the influence of interstitial fluid by
incorporating: (1) strain rate effects, where viscosity effects interact
with the deforming solid skeleton, influencing local mechanical
stimuli; (2) pressure-dependent deformations, allowing for an indirect
representation of fluid-induced mechanical signals; (3) pore evolution,
which modifies the permeability and mechanical response of the system
over time. This approach allows us to evaluate how changes in porosity
and loading conditions influence both the spatial and temporal
evolution of mechanical signals in the system.}

{Previous models for bone remodelling have evolved to incorporate
mechanical stimuli, poroelastic effects, and scaffold bioresorption,
progressively improving their ability to simulate real physiological
conditions. This study builds on two key prior works: one that
introduced a 2D bone/graft model with a diffusive stimulus approach for
bone remodelling but did not extend to full 3D effects or dynamic
loading~\cite{Allena2023}, and another one that analyzed the influence
of loading frequency on bone remodelling within a 3D bone structure,
but without incorporating a bioresorbable scaffold~\cite{Allena2024}.
Although these studies provided valuable information, they did not
address the combined effects of cyclic mechanical loading and scaffold
integration in a 3D bone/graft system, which is the focus of the
present work.}

{Indeed, a key novelty of this study lies in the first integration of a
three-dimensional bioresorbable scaffold (BGB) within a trabecular bone
matrix under dynamic loading conditions. Unlike previous models that
were either limited to two dimensions or considered only the bone
response without scaffold interaction, this approach allows for a more
comprehensive and physiologically relevant simulation of scaffold
remodelling dynamics. The study also extends classical bone remodelling
models by incorporating poroelasticity, which enables the interaction
between solid deformation and interstitial fluid flow---an essential
factor in osteocyte mechanotransduction. This approach enables a more
realistic representation of how oscillatory mechanical forces influence
both bone remodelling and graft resorption over time, providing deeper
insight into the role of mechanical stimulation in optimizing scaffold
integration and bone regeneration. This work also investigates how
dynamic mechanical loads distribute stress and strain within the
bone/graft system, rather than relying on static or simplified loading
conditions as in previous studies. By integrating a 3D scaffold,
dynamic mechanical loading, poroelastic effects, and diffusive
remodelling stimuli, this study aims to represent a significant step
forward in computational modelling of bone healing. The results offer
critical information for optimizing both scaffold design and
rehabilitation protocols, contributing to advancements in orthopedic
applications and tissue engineering strategies.} By incorporating both
bone and graft materials into the model, this study can provide
valuable information about how graft mechanical properties and porosity
influence the bone healing process. The diffusion of the mechanical
stimulus is crucial, as it enables the signal to reach regions far from
the immediate site of the load, promoting a coordinated response
between bone formation and graft resorption. This is particularly
relevant in clinical contexts where bio-resorbable scaffolds are used
to support bone healing after injury or surgery.

The findings of this study improve our understanding of bone/graft
interactions, highlighting the role of mechanical stimuli in driving
the remodelling process in both bone and graft materials. By simulating
a 3D system, the model can capture the complex interaction between
these materials under dynamic mechanical loads, providing a powerful
tool to optimize scaffold design and improve clinical outcomes in bone
repair procedures~\cite{Mullender,cowin1999bm,Zhang}.

\section{The model}\label{The model}

Similarly to our previous works~\cite{Allena2023, Scerrato}, let us
consider a 3D system composed of trabecular bone tissue implanted with
a bio-resorbable scaffold. Both bone and graft materials exhibit
porosity: lacunar-canalicular, inter-trabecular, or collagen-apatites
porosity for the bone and interconnected inner cavities for the graft.
Consequently, we treat them as porous continua and their mechanical
behaviour is governed by the principles of
poroelasticity~\cite{cowin1999bm}
and~\cite{dellIsola2003je,dellIsola1997je,dellIsola2009ijss,burson2023mmcs,madeo2013jmps}. 

The kinematical descriptors that characterize our model are the
displacement $\mathbf{u}$ and the current Lagrangian porosity $\phi$,
which are defined as follows:
{\begin{eqnarray*} 
\mathbf{u} &amp;=&amp; \mathbf{x} - \mathbf{X} \\ 
\phi (\mathbf{X},t) &amp;=&amp; n[\chi(\mathbf{X},t)]J(\mathbf{X},t) 
\end{eqnarray*}}\unskip
where $\mathbf{X}$ is any particle of the system in the reference
configuration, $\mathbf{x}=\chi(\mathbf{X},t)$ is the position of the
particle $\mathbf{X}$ in the current configuration,
$n[\chi(\mathbf{X},t)]$ is the Eulerian porosity and 
$J=\det(\mathbf{F})=\det(\nabla \chi(\mathbf{X},t))$~\cite{coussy2004}.

Throughout evolution, the stiffness of the bone/graft system may attain
significantly softer values than the initial measurements, resulting in
a deformed structure that differs from the reference state. To ensure
consistency, we incorporate nonlinear behaviour most straightforwardly
while assuming that the system remains within the elastic response
range. Furthermore, the rationale for this choice is attributed to the
intricacy of bone microstructure, specifically the trabeculae, which
can exhibit non-linear behaviours even within a minor strain regime, as
evidenced in~\cite{morgan2001jbe}.

{The adopted poroelastic model describes the interaction between a
porous elastic solid and an interstitial fluid. It accounts for fluid
flow through the porous structure and the resulting mechanical
deformations. In the context of bone remodelling, this formulation
allows us to take into account permeability effects; the model
naturally incorporates porosity and permeability, which influence fluid
flow in the lacunar-canalicular network. Increased porosity typically
improves permeability, leading to changes in interstitial fluid flow
and the shear stress sensed by osteocytes. Furthermore, such a model
captures how external mechanical loads induce pressure gradients that
drive fluid through the network, generating mechanical stimuli critical
for osteocyte signaling and bone remodelling.}

The mechanical response of the system is described by the finite-strain
tensor $\mathbf{E}_{ij}(\mathbf{X},t)$ 
{\begin{eqnarray*} 
\mathbf{E}_{ij}(\mathbf{X},t) =
\tfrac{1}{2}(u_{i,j}+u_{j,i}+u_{i,k}u_{k,j}) 
\end{eqnarray*}}\unskip
and the deformation measure for porosity
$\zeta(\mathbf{X},t)$~\cite{Biot1941} reads
{\begin{eqnarray*} 
\zeta(\mathbf{X},t)=\phi(\mathbf{X},t)-\phi^*(\mathbf{X},t)
\end{eqnarray*}}\unskip
with $\phi(\mathbf{X},t)$ and $\phi^*(\mathbf{X},t)$ the Lagrangian
porosities in the current and reference configurations, respectively. 
{Therefore, $\zeta$ can be associated with the change in the fluid
content between the reference to the current configuration.}

Although the system consists of two components---the bioresorbable
graft and the bone---it can be examined as a whole using mixture
theory. To be more precise, the bone exhibits a solid phase as well as
a porous region that is usually occupied by interstitial fluids, blood,
bone cells, and bone marrow. The graft exhibits porosity as well; upon
implantation, its pores are immediately filled with living
things~\cite{lekszycki2012zm,lu2017mms,bednarczyk2016zamp}. Despite
having relatively identical starting conditions, the two portions grow
differently over time. In actuality, three phases can coexist in the
graft during the remodelling process: bioresorbable material, bone, and
physiological fluid.

Accordingly, $\phi^*(\mathbf{X},t)$ can be written as
{\begin{equation}
\phi^*(\mathbf{X},t) = 1 - 
\left(\dfrac{\rho^*_{\mathrm{b}}(\mathbf{X},t)}{\hat{\rho}_{\mathrm{b}}} + 
\dfrac{\rho^*_{\mathrm{g}}(\mathbf{X},t)}{\hat{\rho}_{\mathrm{g}}}\right)
\end{equation}}\unskip
where ${\rho^*_{\mathrm{b}}(\mathbf{X},t)}/{\hat{\rho}_{\mathrm{b}}}$ and
${\rho^*_{\mathrm{g}}(\mathbf{X},t)}/{\hat{\rho}_{\mathrm{g}}}$ are the volume fractions
of the constituents and $\rho^*_{\mathrm{b}}$ and $\rho^*_{\mathrm{g}}$ are the apparent mass
densities of bone tissue and graft artificial material in the reference
configuration, respectively. The superimposed hat indicates the mass
densities evaluated considering the precise volume occupied by those
phases. Since the apparent mass densities of different phases depend on
environment mechanical excitations and time, we assume that the
reference configuration is linked to the apparent mass density in the
reference state, i.e., when zero stress occurs. We assume that the
reference configuration is related to the apparent mass density in the
reference state, or when zero stress occurs, since the apparent mass
densities of different phases rely on time and mechanical excitations
in the environment.  In consideration of the aforementioned
definitions, the current porosity $\phi$ can be expressed as
{\begin{equation}
\phi = \phi^* + \zeta.
\end{equation}}\unskip

The total energy density of the system is given by
\cite{biot1962ap,cowin1983je}
{\begin{eqnarray}
\mathscr{E} = U_s + U_f + U_{fs} &amp;=&amp;
\underbrace{\tfrac{1}{2} \lambda(\rho^*_{\mathrm{b}},\rho^*_{\mathrm{g}})E_{ii}E_{jj}+\mu(\rho^*_{\mathrm{b}},\rho^*_{\mathrm{g}})E_{ij}E_{jj}}_{U_s} 
\nonumber \\&amp;&amp; 
+ \,
\underbrace{\tfrac{1}{2}K_1(\rho^*_{\mathrm{b}},\rho^*_{\mathrm{g}})\zeta^2}_{U_f} + 
\underbrace{\tfrac{1}{2}K_2(\rho^*_{\mathrm{b}},\rho^*_{\mathrm{g}})\zeta E_{ii}}_{U_{fs}}
\end{eqnarray}}\unskip
where all material parameters are considered as a function of the
apparent mass density of the bone and the graft evaluated in the
reference configuration. The subscripts represent the solid part ($s$),
the fluid part ($f$), and the fluid--solid coupling ($fs$) in the given
order. Such an expression is the simplest when considering bone and
graft as isotropic materials (see, e.g.,~\cite{olive2019mmcs} for
generalization). 

$\lambda$ and $\mu$ are the Lam\'{e} coefficients and read
{\begin{equation}
\lambda = \dfrac{\nu Y(\rho^*_{\mathrm{b}},\rho^*_{\mathrm{g}})}{(1+\nu)(1-2 \nu)}, \quad
\mu = \dfrac{Y(\rho^*_{\mathrm{b}},\rho^*_{\mathrm{g}})}{2(1+\nu)}
\end{equation}}\unskip
with the overall Young modulus:
{\begin{equation}
Y = Y_{\mathrm{b}}^{\max} \left(\dfrac{\rho^*_{\mathrm{b}}}{\hat{\rho}_{\mathrm{b}}} \right )^2 + 
Y_{\mathrm{g}}^{\max} \left(\dfrac{\rho^*_{\mathrm{g}}}{\hat{\rho}_{\mathrm{g}}} \right)^2
\end{equation}}\unskip
where $Y_{\mathrm{b}}^{\max}$ and $Y_{\mathrm{g}}^{\max}$ are the maximal values of the Young
moduli for the bone and the graft, respectively. For simplicity, the
Poisson ratio can be assumed to be constant and given. Conversely, we
should consider two evolution rules for the Young moduli and the
Poisson ratios, as discussed in~\cite{branecka2022b,branecka2024cmat}
or, in a broader context, in~\cite{giorgio2023bmm}. This approach
enhances accuracy, but also increases complexity and, in many contexts,
this level of detail is not even necessary. $K_1$ is a compressibility
coefficient defined as the volume of fluid released from unit bulk
volume per unit decrease in pore pressure under the condition of
constant confining stresses~\cite{biot1962ap, Scerrato}. $K_2$ couples
the micro-structure due to pores and the bulk
solid~\cite{Scerrato,biot1962ap}. The expressions of $K_1$ and $K_2$
can be found in~\cite{Scerrato}. As mentioned in~\cite{Scerrato,
Allena2023}, the fluid that fills the pores, the solid matrix of the
bone, the contact between the solid and fluid phases, and the interface
between the graft and the bone are some possible locations for
dissipation in the system. All these sources are taken into account via
a Rayleigh functional as follows:
{\begin{equation}
2\,\mathscr{D}_s = 2\,\mu^{v} (\dot{E}_{ij}\dot{E}_{ij} - 
\tfrac{1}{3}\dot{E}_{ii}\dot{E}_{jj}) + k^{v}\dot{E}_{ii}\dot{E}_{jj}
\end{equation}}\unskip
employing a Kelvin--Voigt model, where $\dot{E}$ is the solid--fluid
rate of deformation and $\mu^{v}$ and $k^{v}$ are two viscous
coefficients.

{Bone remodelling occurs over multiple time scales. Two-time scales are
commonly involved in the overall process. The first is a fast-varying
time scale associated with mechanical loading cycles, such as those
generated by walking or exercising. These loads typically oscillate
with frequencies in the range of seconds to
minutes~\cite{Heinemann2017pe} and induce instantaneous mechanical
responses in the bone/graft system. The second is a slow-varying time
scale that governs the long-term adaptation of bone density and
porosity due to mechanotransduction processes, which evolve over weeks
to months~\cite{eriksen2010remd}. In this paper, we adopt an approach
that follows a well-separated time scale assumption, where rapid
oscillations of the mechanical load contribute to the mechanical
stimulus indirectly affecting the evolution of the long-term parameters
of the material. This allows us to analyze bone remodelling dynamics
without resolving each individual loading cycle, consistent with
previous studies that employ similar temporal separation strategies in
mechanobiological modelling.} Therefore, only the slow-varying-time
scale is considered, as the rapid oscillations due to the former
fast-varying-time scale can be treated from the perspective of
well-separated time scales, as explained in~\cite{scerrato2022s}. 

The Generalized Principle of Virtual Work is as follows:
{\begin{equation}
\int _{\mathscr{B}^{*}} \delta\mathscr{E}\,\mathrm{d}\mathscr{B}^{*} + 
\int _{\mathscr{B}^{*}} \dfrac{\partial\mathscr{D}_{s}}{\partial\dot{E}_{ij}}\delta E_{ij} \partial \mathscr{B}^{*} = 
\int _{\mathscr{B}^{*}} \delta\mathscr{W}^{\mathrm{ext}}\,\mathrm{d}\mathscr{B}^{*}
\end{equation}}\unskip
where $\delta\mathscr{W}^{\mathrm{ext}}$ is the virtual work generated
by the external loads, which can be written as
{\begin{equation}
\delta \mathscr{W}^{\mathrm{ext}} = 
\int _{\partial_{\tau} \mathscr{B}^{*}}\tau_{i} \delta u_{i} \,\mathrm{d} \mathscr{S}^{*} +
\int _{\partial_{\Xi} \mathscr{B}^{*}} \Xi \delta \zeta \,\mathrm{d}\mathscr{S}^{*}.
\end{equation}}\unskip
In this equation the first term indicates the forces per unit surface
(i.e., $\tau_{i}$ on the part of the boundary $\delta_{\tau}
\mathscr{B}^{*}$) and the second term relates to the effects of a pore
pressure $\Xi$ on the boundary $\delta_{\Xi} \mathscr{B}^{*}$. In this
study, since we focus solely on the slow-varying time scale, we can
neglect the effects of inertia. To fully understand these effects (as
shown in similar complex structures, see,
e.g.,~\cite{eugster2022mmcs,placidi2024mrc}), it would be necessary to
introduce a fast-varying time scale. However, as previously remarked,
for the sake of simplicity, we have chosen not to include that in this
work. 

Once the graft has been implanted, bone remodelling takes place around
it and, after a while, also within it. We describe the evolution over
time of the apparent mass density of the bone ($\rho_{\mathrm{b}}$) and the
graft ($\rho_{\mathrm{g}}$) with respect to the reference configuration as
follows:
{\begin{equation}
\begin{cases}
\dfrac{\partial \rho_{\mathrm{b}}^{*}}{\partial  t}(\mathbf{X},t) = 
A_{\mathrm{b}}[S(\mathbf{X},t)]H[\phi(\mathbf{X},t)] &amp; \mbox{with }
0&lt;\rho_{\mathrm{b}}^{*}\leq\hat{\rho}_{\mathrm{b}} \vspace*{6pt}\\
\dfrac{\partial \rho_{\mathrm{g}}^{*}}{\partial t}(\mathbf{X},t) = 
A_{\mathrm{g}}[S(\mathbf{X},t)]H[\phi(\mathbf{X},t)]  &amp;\mbox{with }
0&lt;\rho_{\mathrm{g}}^{*}\leq\rho_{\mathrm{g}}^{*}(\mathbf{X},0)
\end{cases}
\end{equation}}\unskip
where $S(\mathbf{X},t)$ is the mechanical stimulus and $A_{\mathrm{b}}$ and $A_{\mathrm{g}}$
are two functions assumed to be piecewise linear as follows: 
{\begin{equation}
A_{\{\mathrm{b},\mathrm{g}\}}[S(\mathbf{X},t)] = \begin{cases}
s_{\{\mathrm{b},\mathrm{g}\}}[S(\mathbf{X},t) - S_s] &amp; \mbox{for } (S - S_s)\geq 0\\
0\ &amp; \mbox{for } S_r &gt; S &gt; S_s \\
r_{\{\mathrm{b},\mathrm{g}\}} [S(\mathbf{X},t)-S_r] &amp; \mbox{for } (S - S_r)\leq 0 .
\end{cases}
\end{equation}}\unskip
{The notation $\cdot_{\{\mathrm{b},\mathrm{g}\}}$ has been introduced to combine two
formulas into a single, unified expression. The coefficients $s$ and
$r$ represent the rates of growth and resorption, respectively, with
the subscripts ${\mathrm{b},\mathrm{g}}$ denoting bone and graft, respectively. The
quantities $S_s$ and $S_r$ serve as two thresholds that define the
``lazy zone'', which corresponds to the homeostatic state. These
reference values of the stimulus are assessed with the understanding
that they are energy-based. Consequently, we establish an assumed
typical range of strain between 600 and 2000 microstrains from which we
subsequently calculate the corresponding energy associated with that
level of deformation.}  {It has to be noticed that $s_{\mathrm{g}}=0$. In fact,
we assume that the scaffold material degrades over time without
undergoing new formation, whereas bone can grow in response to
mechanical stimuli ($s_{\mathrm{b}}&gt;0$). This assumption is consistent with the
behavior of most bioresorbable materials, which are designed to be
gradually resorbed and replaced by new bone tissue rather than
self-renewing. Additionally, we have set $r_{\mathrm{g}}=r_{\mathrm{b}}$. We are aware
that kinetics of bone resorption and scaffold degradation are governed
by distinct biological and physicochemical processes. Bone resorption
is primarily regulated by osteoclast activity, whereas scaffold
resorption depends on several factors (i.e.\ hydrolysis, enzymatic
degradation\,\ldots). Nonetheless, this assumption ensures that resorption
occurs at comparable rates in both materials, facilitating a direct
comparison of their evolution under mechanical stimulation. In
addition, it is a good practice that when designing the graft, the
properties of the inserted material closely align with those of the
original bone tissue. This adherence is crucial to ensure the optimal
functionality of the replacement material. A significant disparity in
the behavior of the graft compared to the native bone can lead to
complications in both the short- and long-term reconstruction
processes.}

In this work, we adopt a classical approach grounded in
phenomenological evolution rules. However, this simplification can be
easily improved through variational methods, as outlined
in~\cite{grillo2023mmcsa,grillo2023mmcsb}.

During the remodelling process, porosity is crucial and is taken into
account by the function $H[\phi(\mathbf{X},t)]$. The pores facilitate
the access of bone cells to the more remote regions of the system,
allowing them to resorb or generate bone tissue in reaction to external
forces. Thus, two potential results arise: (i) inadequate porosity
hinders remodelling because no room is left for the bone cells to
operate, or (ii) excessive porosity fails to support bone cells, both
resulting in unchanging densities. Consequently, we have assumed that
porosity reaches a maximum value of $H$ at approximately 0.6, which
provides optimal conditions for bone remodelling when porosity is
distanced from the two extreme scenarios~\cite{giorgio2019bmm}. 

Traditionally, osteocytes perceive the mechanical stimulus
$S(\mathbf{X},t)$, allowing the bone to assess and modify its overall
state. 

Here, we assume that the solid phases---and, by extension, the
mechanical properties---are altered by both osteoclasts and
osteoblasts. As we have previously proposed
in~\cite{giorgio2019bmm,Allena2023,Allena2024}, the stimulus is defined
via a diffusion equation. 
{\begin{equation}
\frac{\partial S}{\partial t} - \kappa \Delta S +\beta S =f(\rho_{\mathrm{b}}^*) 
\mathscr{E}
\end{equation}}\unskip
with $\kappa$ and $\beta$ being two scalars;  {$\beta$ is only
activated when the stimulus is produced ($S&gt;0$), and it is associated
with a metabolic resorption of the stimulus that prevents excessive
accumulation of the stimulus.}   {Here, we assume that $\kappa$ and
$\beta$ are the same in scaffold and bone. Nonetheless, we are aware
that the actual properties of these two materials may vary in terms of
porosity and microstructure. This simplification allows us to focus on
the basic phenomena of dynamic loading and mechanotransduction.} The
function $f(\rho_{\mathrm{b}}^*)$ is a proper signalling weight as defined
in~\cite{Allena2024} and reads
{\begin{equation}
f(\rho^*) = \dfrac{1}{2}\left \{\tanh\left[ 20\left
(\dfrac{\rho^*}{\hat{\rho}} -0.15 \right) \right]
+\tanh [20(1-0.15)]\right\}.
\end{equation}}\unskip

This function depends on the apparent mass density of the tissue and
increases monotonically. Assuming, as a first approximation, that
osteocytes are distributed nearly equally throughout the tissue, we can
infer that their density per unit volume is proportionate to the volume
fraction.

The stimulus is equal to $0$ and 367.57~J in the graft and in the
bone, respectively, at the initial time $t=0$.  {In the graft, the
initial stimulus $S$ is zero because at the beginning there are no
osteocytes capable of producing it. In bone tissue, the initial
stimulus is set within the range of the homestatic state ($S_r &gt; S &gt;
S_s$); therefore, we start from an equilibrium state.} Additionally, no
stimulus is exchanged with the outside along the boundaries of the
system. 

In general, in the biological system considered, the accuracy of the
models is heavily based on the precise estimation of parameters, which
represent rates, thresholds, or interactions within the system. Given
the inherent complexity and the potential lack of direct measurability,
parametric evaluation through proper identification procedures becomes
essential. These procedures, often involving optimization techniques,
statistical inference, or other ad hoc methods, aim to reconcile models
with experimental or observational data. By systematically refining
parameter estimates, such methods improve model reliability and
predictive power, allowing deeper insight into the underlying
biological processes and better handling of uncertainties intrinsic to
these
systems~\cite{valmalle2023ejma,aydin2022,yang2022,fedele2024cmt,misra2021zamp}. 
{In this study, we did not perform a formal parameter identification
procedure. Instead, the model parameters were estimated through
numerical calibration, ensuring that the simulated response captured
key qualitative aspects observed in experimental studies. Our primary
objective at this stage was to investigate the qualitative influence of
mechanical loading on bone remodelling, rather than to achieve a
precise quantitative fit with experimental data. This approach allowed
us to explore fundamental trends, assess the sensitivity of the model
to different mechanical stimuli, and identify relevant parameter ranges
for future refinement. A more rigorous parameter identification
process, which incorporates the fitting of experimental data and the
quantification of uncertainties, is planned for subsequent studies to
enhance the predictive accuracy of the model.}

\section{Numerical implementation for an illustrative case}
\label{Numerical implementation}

In this work, the numerical investigations were carried out using the
COMSOL Multiphysics software. A weak formulation that applies finite
element analysis was used to directly solve the partial differential
equations that constitute the previously described mathematical model.
We employ this direct approach to create a sufficiently accurate model
that utilizes the general-purpose capabilities of the chosen software,
even if it may not be highly efficient. In the second stage of
development, we can apply numerical procedures that are appropriately
optimized for the model generated in the previous steps (see,
e.g.,~\cite{greco2021cmame,tran2020mmcs}).

We consider a cylinder of length $L$ and radius $r$, which represents a
typical femoral diaphysis. The specimen is composed of bone at the
bottom and at the top and by a graft of length $L_{\mathrm{g}}$ in the middle,
resulting in a BGB structure (Figure~\ref{Fig1}). The cylinder is fully
constrained at the bottom, while on the top surface, a sinusoidal force
producing torsion and a pressure are applied as follows:
{\begin{equation}
f = (X_2\mathbf{e}_{\boldsymbol{2}} - X_1\mathbf{e}_{\boldsymbol{1}}) F_0 
\sin(2 \uppi \Omega t) - [F_1 + F_2  \sin (2 \uppi \Omega t)]
\mathbf{e}_{\boldsymbol{3}} 
\end{equation}}\unskip
where $F_0$, $F_1$ and $F_1$ are nominal forces, $\Omega$ is the
frequency, and
$\mathbf{e}_{\boldsymbol{1}},\mathbf{e}_{\boldsymbol{2}}$ and
$\mathbf{e}_{\boldsymbol{3}}$ are the unit vectors of the coordinate
system.  {We remark that $\Omega$ denotes a frequency that gradually
changes over a period of weeks to months and does not capture the rapid
frequency at the mechanical load time scale, consistent with the
assumption of well-separated time scales. In fact, while the applied
load in our model is oscillatory (at the slow-varying time scale), its
influence on tissue adaptation is considered through the cumulative
mechanical stimulus, which integrates the effects of strain energy
density and fluid--solid interactions over time. In other words, under
our assumption, frequency load pertains exclusively to the slow time
scale and thus can be interpreted as the cyclical process of loading
and unloading over an extended period. Within the context of the
well-separated time scale assumption, the frequency at the fast time
scale is not considered, as it only indirectly influences the system
through its effect on the effective amplitude of the mechanical load
(see for more details~\cite{scerrato2022s}).} The model parameters are
reported in Table~\ref{table}.

\begin{figure}
\includegraphics{fig01}
\caption{\label{Fig1}Schematic representation of the 3D system under
study. The sample is constituted by two bone (B) regions and a graft
(G) region in the middle in a BGB sequence. The sample is fully clamped
at the bottom, while sinusoidal torque and pressure are applied to the
top surface.}
\end{figure}

%tab1
\begin{table}
\caption{\label{table}Structural parameters of the model used in the
numerical simulations}
\begin{tabular}{cc}
\tbody
$L$ (cm) &amp; 15\\
$r$ (cm) &amp; 2\\
$L_{\mathrm{b}}$ (cm) &amp; 5\\
$L_{\mathrm{g}}$ (cm) &amp; 5\\
$\hat{\rho}_{\mathrm{b}}$ (kg/m$^{3}$) &amp; 1800\\
$\hat{\rho}_{\mathrm{g}}$ (kg/m$^{3}$) &amp; 1800\\
$\nu$ &amp; 0.3\\
$Y_{\mathrm{b}}^{\max}$ (GPa) &amp; 18\\ 
$Y_{\mathrm{g}}^{\max}$ (GPa) &amp; 14.4\\ 
$\nu^{v}$ (N$\cdot$s/m$^{2}$) &amp; $1 \times 10^{-4}$\\
$k^{v}$ (N$\cdot$s/m$^{2}$)&amp; $1 \times 10^{-4}$\\
$s_{\mathrm{b}}$ (s/m$^{2}$) &amp; $4.9 \times 10^{-9}$\\
$r_{\mathrm{b}}$ (s/m$^{2}$) &amp; $6.13 \times 10^{-9}$\\
$r_{\mathrm{g}}$ (s/m$^{2}$) &amp; $6.13 \times 10^{-9}$\\
$\kappa$ (m$^{2}$/s)&amp; $1.6 \times 10^{-4}$\\
$\beta$ &amp; 0.1\\
$\Omega$ (CpT) &amp; 0.01, 0.1, 1\\
$F_0$ (MPa/m) &amp; 200 \\
$F_1$ (MPa) &amp; 0.16 
\botline
\end{tabular}
\end{table}
 
The numerical simulations have been performed through a finite element
analysis using COMSOL Multiphysics.

\section{Results and discussion}\label{Results and discussion}

To analyze our numerical results, we consider six points inside the
cylindrical specimen. Three of them ($P_1, P_2, P_3$) are located along
the vertical axis of the cylinder at a height of 2.5, 7.5 and 12.5~cm.
The points $P_4$, $P_5$, and $P_6$ are moved by 1.5~cm along the
$\mathbf{e}_2$ axis to have probe points near the lateral surface of
the sample (Figure~\ref{Fig1}). 

In the first series of simulations, we have evaluated the impact of the
frequency $\Omega$ on the system when $F_0 = 200$~MPa/m,  $F_1 =
0.16$~MPa and $F_2 = 0.016$~MPa. {These force amplitudes are
numerically calibrated to produce a level of stimulus that stimulates
the generation of new bone tissue. In essence, we consider the scenario
where the stimulus $S$ exceeds $S_s$, the synthesis threshold.
Specifically, the maximum torque applied has an amplitude of 50.26~Nm,
which is compatible with a reference value taken from the literature
(see, e.g.,~\cite{kotzar1995jor}).} 

The first set of graphs in Figure~\ref{Fig3} presents the evolution of
the stimulus $S$ over time at three different frequencies, specifically
$\Omega = 0.01$ Cycles per unit Time (CpT) (Figure~\ref{Fig3}a,d),
$\Omega = 0.1$ CpT (Figure~\ref{Fig3}b,e) and $\Omega = 1$ CpT
(Figure~\ref{Fig3}c,f). At low frequency (i.e., $\Omega = 0.01$ CpT,
Figure~\ref{Fig3}a,c), the stimulus response across all points is
relatively flat, showing minimal growth, particularly at points $P_1,
P_3$ and $P_4, P_6$, represented by the red lines. These points
experience very low stimuli, suggesting either less mechanical
interaction or areas of low strain, potentially in regions further from
the loading centre or areas with less bone density. At intermediate
frequency (i.e.~$\Omega = 0.1$~CpT, Figure~\ref{Fig3}b,e), the stimulus
increases more significantly, with a sharp rise over time at certain
points (notably $P_1, P_4$ and $P_3, P_5$, marked by blue and red
lines, respectively). This response suggests that intermediate
frequencies may optimize the mechanical environment to trigger a more
substantial biological response. The green lines referring to the
points in the graft (points $P_2$ and $P_5$) also show mild growth,
indicating that these points are moderately responsive to frequency
changes, though less active than the blue lines. This is because the
stimulus originates from the osteocytes in bone tissue and then also
reaches the graft with a slight delay. Finally, high-frequency loading
(i.e., $\Omega = 1$ CpT, Figure~\ref{Fig3}c,f) leads to a much sharper
increase in stimulus at points $P_1, P_3$ and $P_4, P_6$ (blue and red
lines), suggesting that under these conditions, those points respond
more significantly to mechanical stimulus. This could reflect a shift
in the mechanical load distribution at higher frequencies, where, even
initially, the less active points become more stimulated.

\begin{figure}
{\vspace*{-6pt}}
\includegraphics{fig02}
{\vspace*{-6pt}}
\caption{\label{Fig3}Influence of the frequency $\Omega$ on the
stimulus $S$ for points $P_1, P_2, P_3$ (a,b,c) and points $P_4, P_5,
P_6$ (d,e,f). (a,d) $\Omega = 0.01 $, (b,e) $\Omega = 0.1$ CpT, (c,f)
$\Omega = 1$ CpT when $F_0 = 200$~MPa/m and $F_1 = 0.16$~MPa.}
{\vspace*{-6pt}}
\end{figure}

Figure~\ref{Fig4} explores the evolution of the {normalized} bone mass
density $\rho_{\mathrm{b}}^*$ over time, representing the gradual adaptation of
bone tissue in response to the mechanical stimuli seen in
Figure~\ref{Fig3}. The same points are analyzed under the three
frequency regimes. Figure~\ref{Fig4} panels (a) and (d) demonstrate a
minimal change in bone mass density, with the points $P_1, P_3$ and
$P_4, P_6$ (blue and red lines, respectively) maintaining the initial
density value (i.e., $\rho_{\mathrm{b}}^* = 0.5$), while bone density at the
points $P_2$ and $P_5$ in the graft remains nearly equal to zero. This
indicates that, at low frequencies, remodelling is minimal and bone
mass density remains relatively static. At $\Omega = 0.1$ CpT
(Figure~\ref{Fig4}b,e), bone mass density increases markedly over time
for points $P_1, P_3$ and $P_4, P_6$ (blue and red lines,
respectively). The points $P_2$ and $P_5$ (green lines) show little or
no change, suggesting that these regions might experience less strain
or be further from the main loading sites. Similarly to intermediate
frequencies, Figure~\ref{Fig4}c and f show a pronounced increase in
bone density for points $P_1, P_3$ and $P_4, P_6$ (blue and red lines,
respectively) at $\Omega = 1$ CpT. This behaviour implies that
high-frequency loading could significantly enhance the remodelling
process, resulting in higher bone mass density. The nearly flat green
lines (points $P_2$ and $P_5$) suggest that the regions associated with
these points are less responsive even under high-loading conditions.
Since this \mbox{region} is associated with the graft, this means that the
production of the stimulus that is relevant in bone tissue is not yet
sufficient to colonize the implant region and therefore has a relevant
effect. Such results can also be observed in Figure~\ref{Fig2}a,b and
c, where we show the distribution of the  {normalized} bone mass
density $\rho_{\mathrm{b}}^*$ in the 3D BGB structure at $t = 25$ weeks. One can
observe that even with relatively low forces, increasing the frequency
significantly enhances bone formation. This indicates that the cyclic
nature of stimulation plays a crucial role in triggering remodelling
responses in bone tissue.

\begin{figure}
\includegraphics{fig03}
\caption{\label{Fig4}Influence of the frequency $\Omega$ on the 
{normalized} bone mass density $\rho_{\mathrm{b}}^*$ for points $P_1, P_2, P_3$
(a,b,c) and points $P_4, P_5, P_6$ (d,e,f). (a,d) $\Omega = 0.01 $,
(b,e) $\Omega = 0.1$ CpT, (c,f) $\Omega = 1$ CpT when $F_0 = 200$~MPa/m
and $F_1 = 0.16$~MPa.}
\end{figure}

\begin{figure}
\includegraphics{fig04}
\caption{\label{Fig2}Normalized bone density $\rho_{\mathrm{b}}^2$ at $t = 25$
weeks for $F_0 = 200$~MPa/m, $F_1 = 0.16$~MPa and $\Omega = 0.01$ CpT
(a), $\Omega = 0.1$ CpT (b), $\Omega = 1$ CpT (c) and $F_0 =
1000$~MPa/m, $F_1 = 0.16$~MPa and $\Omega = 0.1$ CpT (d).}
\end{figure}

Figure~\ref{Fig5} focuses on the evolution of the {normalized} density
of the graft mass $\rho_{\mathrm{g}}^*$ at the same points at different
frequencies. Unlike bone mass density, the density of the graft mass
appears to decrease across all points regardless of frequency. For low,
intermediate and high frequencies, Figure~\ref{Fig5} panels from (a) to
(f) show a sharp decline in the density of the graft mass within the
first few weeks. The green lines (points $P_2$ and $P_5$) fall to near
zero after approximately 10--15~weeks. This behaviour suggests a process
where the graft resorbs over time, potentially being replaced by newly
formed bone tissue in a process characteristic of osteointegration, as
it can be observed at intermediate and high frequency 
(Figure~\ref{Fig4}b,c,e,f). This type of study is crucial for assessing
the timing of resorption and subsequent bone formation. In fact, if the
graft is completely resorbed prior to bone tissue formation, the
mechanical integrity of the system could be severely compromised.

\begin{figure}
\includegraphics{fig05}
\caption{\label{Fig5}Influence of the frequency $\Omega$ on the graft
{normalized} mass density $\rho_{\mathrm{g}}^*$ for points $P_1, P_2, P_3$ (a,b,c)
and points $P_4, P_5, P_6$ (d,e,f). (a,d) $\Omega = 0.01$ CpT, (b,e)
$\Omega = 0.1$~CpT, (c,f)~$\Omega = 1$~CpT when $F_0 = 200$~MPa/m and
$F_1 = 0.16$~MPa.}
\end{figure}

To enhance osteointegration, we have increased the mechanical loads.
While we observed that the pressure does not play any role even at
higher values, torsion seems to affect the overall response of the BGB
structure. Here, we report the results of a second series of
simulations where $F_0 = 1000$~MPa/m, $F_1 = 0.16$~MPa and $\Omega =
0.1$~CpT.

Figure~\ref{Fig6} shows the temporal evolution of the stimulus $S$, the
bone mass density $\rho_{\mathrm{b}}^*$, and the graft mass density $\rho_{\mathrm{g}}^*$ for
specific points $P_1,P_2,P_3$ (in plots a,b,c) and $P_4,P_5,P_6$ (in
plots d,e,f). The curves in Figure~\ref{Fig6}a and d exhibit a
smooth increase of the stimulus $S$ for the points $P_1$, $P_3$ and
$P_4$, $P_6$ (blue and red curves, respectively), while in the graft
(points $P_2$ and $P_5$, green lines) $S$ undergoes a steep upward
trend after 10 weeks. Such values correlate with increased bone
remodelling activity, which is consistent with what we see in bone
density plots (Figure~\ref{Fig6}b,e). On the one hand, both in
the bone and in the graft, $\rho_{\mathrm{b}}^*$ increases, reaching values of
about 0.95 kg/m$^3$ and 0.85 kg/m$^3$ for points $P_1$, $P_3$, $P_4$,
$P_6$ (blue and red lines) and $P_2$, $P_5$ (green lines),
respectively. On the other hand, the density of the mass of the graft
$\rho_{\mathrm{g}}^*$ is reduced at all points, meaning that the graft material is
gradually resorbed or remodelled as bone tissue develops. Furthermore,
depending on the level of the stimulus and its distribution over the
sample, the mass density of the graft may remain different from zero
before the generation of new bone tissue in the same region to prevent
a possible rupture.

\begin{figure}
{\vspace*{-5pt}}
\includegraphics{fig06}
{\vspace*{-5pt}}
\caption{\label{Fig6}Evolution over time of the stimulus $S$, the bone
{normalized} mass density $\rho_{\mathrm{b}}^*$ and the graft  {normalized} mass
density $\rho_{\mathrm{g}}^*$ for the points $P_1, P_2, P_3$ (a,b,c) and $P_4,
P_5, P_6$ (d,e,f) when $F_0 = 1000$~MPa/m, $F_1 = 0.16$~MPa and $\Omega
= 0.1$ CpT.}
{\vspace*{-5pt}}
\end{figure}

Overall, one can conclude that low-force and low-frequency loading
result in minimal changes in bone density. At the same time,
high-frequency and high-force conditions promote optimal bone growth
and graft integration. Therefore, mechanical stimuli must be carefully
calibrated in terms of both force and frequency to maximize bone growth
and graft integration. For practical applications in bone grafts or
orthopaedic implants, these results could guide the selection of
rehabilitation protocols or the design of mechanical devices that apply
controlled cyclic loading to encourage efficient bone regeneration and
implant integration.

\section{Conclusion}

This study presents a 3D poromechanical model for simulating bone
remodelling in systems incorporating bio-resorbable grafts, providing a
detailed analysis of the interplay between mechanical loads and
remodelling dynamics.  {Our results demonstrate that the application of
dynamic mechanical loads in a 3D BGB structure significantly influences
both bone formation and graft resorption, emphasizing the importance of
mechanical stimulation in scaffold integration. Unlike previous models
that considered either dynamic loading in simplified geometries or
studied 3D remodelling under static conditions, our approach combines
these elements to capture the full spatial and temporal complexity of
bone adaptation. The findings suggest that optimizing load frequency
and intensity in a 3D environment can lead to improved bone
regeneration, which has direct implications for the design of
bioresorbable scaffolds and post-surgical rehabilitation strategies.}
More specifically, the simulations reveal that higher frequencies and
load magnitudes enhance both bone regeneration and graft integration,
while low mechanical loads yield poor remodelling effects from the
clinical viewpoint.

The findings underscore the importance of calibrating mechanical loads
in terms of intensity and frequency to optimize bone healing and
scaffold integration. This model could serve as a valuable tool for
improving the design of bio-resorbable grafts and developing targeted
rehabilitation protocols. By capturing complex spatial and temporal
interactions in bone/graft systems, the study could contribute to the
advancement of both theoretical understanding and clinical applications
in bone repair and regeneration. In future development, we want to add
the effects of damage~\cite{placidi2022me,espino2023} and
dissipation~\cite{chekeres2024cmt,cuomo2024mmcs} to our model, since we
believe that they play a crucial role in bone remodelling by driving
biological processes that repair microdamage and maintain structural
integrity, ensuring that the bone adapts effectively to mechanical
stresses over time. This aspect will improve our understanding of the
process and will allow us to add a crucial piece to the puzzle of bone
remodelling.

{Furthermore, bone remodelling is governed by cellular processes and
external mechanical loads that occur over the course of minutes, while
the comprehensive effects of this \mbox{mechanobiological} phenomenon become
evident after a duration of several weeks. The simultaneous presence of
these two distinct time scales complicates the design of experiments
aimed at evaluating and quantifying specific mechanical variables.
Therefore, another objective of this paper is to present a modelling
tool with sufficient accuracy for the planning of a measurement
campaign, with the ultimate aim of validation. This methodological
approach is crucial to ensure that experimental data are acquired in a
manner that directly assesses model predictions while minimizing
uncertainties and confounding factors. The use of a model as a guide
can facilitate the validation process by identifying key predictions
from the model, selecting measurable quantities, and optimizing
experimental conditions. This involves determining the physical
quantities predicted by the model (e.g., fluid pressure, strain fields,
porosity effects on flow); identifying the most sensitive parameters,
namely those that significantly influence the outcome and should be
measured with high precision; selecting experimental variables that
directly correspond to model outputs (e.g., pore pressure via pressure
sensors, strain via digital image correlation or strain gauges, fluid
flow via tracer experiments); considering the potential of indirect
measurements (e.g., ultrasound or MRI for fluid flow) to provide the
required data; utilizing the model to simulate various testing
conditions to identify those that maximize the differentiation between
hypotheses; ensuring that the selected conditions (e.g., loading
frequency, amplitude, or boundary conditions) produce measurable
responses; ensuring that experimental data are collected with adequate
spatial and temporal resolution to capture significant model features;
and minimizing measurement redundancy while ensuring ample data for
robust validation.}

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and have
declared no affiliations other than their research organizations.

\section*{Acknowledgments}

The first author is supported by the Institut Universitaire de France.

\back{}

\printbibliography
\refinput{crmeca20240989-reference.tex}

\end{document}
</pre></body></html>