Abstract
Predictive theoretical–numerical modeling of the behavior and evolution of biological tissue is a difficult task since many of the required knowledge tools (experimental, theoretical and numerical) are still not well understood. We present here some methodologies and results specific to multiscale and multiphysics numerical modeling of biological tissues applied to the predictive behavior of cortical veins depending on their local constituents’ microstructure and for bone remodeling and reconstruction as a function of the local mechanobiology. Although further work is required to improve the accuracy of the developed models, the proposed approaches highlight their potential usefulness for understanding the mechanical-biological couplings, short and long term predictions of biological evolutions as well as possible further transfer to medical applications.
Introduction
The use of adequate numerical models to predict the mechanical behavior of soft tissues [1,2] is often critical for patient health and safety. To improve these models up to the medical standards and regulations, more detailed information is still needed [3–6]. In addition, for many medical applications, post-aesthetic [7] and geometry [8] are important factors for patients. For bone reconstruction, numerous questions remain in understanding the specific mechanobiology and nutriment supply [9–12] since it is difficult to quantify the patient biological parameters during and post-surgery. Numerical models can help the surgeon in its everyday job [13–16], but it is also desirable to develop precise 3D models to predict long term evolutions [17–21].
In this work, we describe the development of such numerical models for the prediction of biological tissue evolutions accounting for both the mechanical and biological aspects. Example of the cortical vein microstructure is implemented within the theory of homogenization [1,3,5,22], and using experimental input data [4,6], in order to provide useful information on their potential rupture leading to acute subdural hematoma. In a second stage, we describe the complex mechanobiological interactions [9,23–29] occurring within bone tissue remodeling knowing the biological phenomena in play and following ordinary differential equations slowly varying in time. The actor cells such as osteoclasts and osteoblasts are activated by means of a biological stimulus defined in the form of an energy function based on the applied mechanical loads [16–21,30]. We account for the interplay between biological and mechanical effects known to be important for the bone tissue synthesis and we study the effect of different applied mechanical load intensities and directions to characterize the quality, rate and density of bone remodeling. More specifically, the links existing between the micro and the macro scales are studied. We account for the local cell activities and localization, reaction to local stimulus in order to predict the microstructure evolutions and overall response. This will enhance understanding and optimize the accuracy of the long term predictions for bone implant, remodeling and post-surgery reconstruction.
Modeling soft tissue: The cortical vein behavior
Experiments
The establishment of a precise model of biological soft tissue mechanical behavior requires in depth knowledge of the tissue microstructure since its macroscopic behavior is fully dependent on the mechanical properties of each of its local constituents. However, vein wall microstructure is complex and requires detailed investigation. In the current work, to extract these local constituents’ properties, different scales were identified as presented in Fig. 1. This was done through optical and numerical microscopy and enabled first to localize the cortical vein organization, collagen fibers orientation and wall thickness, and a first evaluation of the vein wall microstructure.

Identification of the different working scales for the experiments and numerical model.
Regarding the mechanical behavior of the complex network of collagen fibers, poor information is currently available in the literature [31,32]. In this work, we used biphotonic confocal microscopy (see Fig. 2(a)) to identify the collagen fiber behavior using a micro tensile test machine placed under the confocal microscope.

Biphotonic confocal microscopy experiment (a) with elongated collagen fibers (b).
Observations were made on a number of specimens onto which mechanical tension was applied. This experiment enabled the observation of the collagen fibers reorientation as a function of the applied deformation ϵ (see Fig. 2(b)) but also the characterization of the associated mechanical behavior through the extraction of the relation force applied vs measured displacement. This evolution of the collagen network also provided useful information for the construction of the multiscale model (see results section).
Modeling the mechanical behavior of cortical veins was done in two steps. The first one consists in representing the localization of each of the cortical vein along the superior sagittal sinus (SSS) as shown on Fig. 3(a) by the definition of an adequate representative volume element (RVE, Fig. 3(b)) and provide an adequate homogenization model along the SSS to optimize the corresponding macroscopic mechanical behavior. Homogenization is done through accounting of the variations of the mechanical parameters of this RVE (viscosity μ, stress σ, deformation ϵ,

Definition of the homogenization model for the cortical vein behavior (reprinted with permission from Mathematical Sciences Publishers, JOMMS 7:6 (2012), 597 – Mathematical Sciences Publishers).
The displacement field, using asymptotic homogenization technique along the SSS is provided by
Where H is a tensor obtained by solving the equilibrium equation. This stiffness tensor
Once the macroscopic behavior identified, we need to integrate the vein local constituents’ mechanical behavior to link them to the macroscopic response of the structure. The second step of the analysis consists in defining the specific microstructure of the vein wall (see Fig. 4) in order to be able to develop the adequate corresponding theoretical model.

Definition of the cortical vein wall microstructural scales.
Accounting for the local constituents’ behavior in the global macroscopic response requires again the use of homogenization technique. Assuming that each of the local constituents have linear elastic behavior following Hooke’s Law, we can write for each constituent i;
Defining the overall mechanical behavior of the cortical vein imposes to couple the macroscopic asymptotic homogenization model with the multiscale microstructure model. Results of each are presented individually here.
Using the asymptotic homogenization model together with local constituents’ individual mechanical properties, we found the following homogenized stiffness coefficients for the SSS-cortical vein junction:
With

Reorganization of the vein wall collagen structure as a function of applied mechanical strain: (left) view through biphotonic confocal microscopy, (right) evolutionary graph of the two collagen fiber families.
Next, the local microscopic effects depending on the vein wall collagen microstructure were identified. This was done both experimentally and theoretically. The collagen fiber distribution was extracted from confocal microscopy experiments (Fig. 2). The determination of the structural reorientation of the collagen fibers as a function of elongation was necessary in order to quantify the corresponding macroscopic mechanical behavior of the vein wall. The ImageJ software [34,35] was used together with the module OrientationJ [36]. Figure 5 (left) shows the visible reorientation of collagen fibers vs applied mechanical strain from 0% to about 70% Henckystrain. We can observe a clear reorientation of the fibers going from globally vertical (on the figure) to nearly horizontal orientation. The evolution on the two main fiber families is shown on Fig. 5 (right). The angular variation goes from an average value of 71.58° down to an average value of 16.87° from the horizontal. The average amplitude of reorientation is about 55°. The reorientation of the collagen fibers provides a precise quantification of the 3D evolution as a function of the applied mechanical force.
We extracted the corresponding mechanical behavior of the structure and implemented it in the developed theoretical numerical model. Finally, the constitutive equation was integrated in the multiscale numerical model to extract the overall vein wall behavior. A comparison of the theoretical multiscale model results with experimental data available in the literature [37] is presented in Fig. 6.

Comparison of numerical model and experiments from [37] for experimental tensile behavior of human parasagittal bridging veins. The stress is the average in the main vein wall direction during tension and Hencky strain is the Cauchy–Green deformation in the tension direction.
A very good correlation between experiment and theoretical numerical model is observed and showed the effectiveness of using these approaches for the prediction of the macroscopic soft tissues deformation.
Context
Bone is a complex biological tissue and the rules governing its evolution are also complex [9,23–29]. Although it is well know that bone density evolve as a function of the applied external mechanical forces [38–40], this evolution is mainly dependent on the bone biology [9,41,42]. Understanding of these phenomena is important for the prevision of bone remodeling and bone reconstruction in many orthopedic and surgical applications. For example, in orthodontics or maxillo facial surgery, the application of mechanical forces has a direct impact on bone remodeling and reconstruction, and will modify the bone morphology accordingly to the local biology (vascularization, bone density, etc
Theoretical numerical model
We recall here the main steps of the developed continuum model. More detailed information can be found in [17–21]. We describe the evolution of a two-solid continuum mixture associating any given material point
The mechanobiological coupling is provided by the biological stimulus S defined on the volume Ω [19] by
Finally, the bone density kinetic is defined through the ordinary differential remodeling equation
This theoretical model was implemented in numerical simulations for different applications and external boundary conditions. We present here the bone density kinetics and microstructures observed on the example cases presented on Fig. 7.

Geometries of numerical models. Left (a): Macroscopic femur head with top applied load, right: (b) microstructure trabeculae with 45° angular load and (c) triangular load.
Case one (Fig. 7(a)) represents a continuous 3D femoral head fixed at its bottom end and loaded mechanically over the top area corresponding to the applied body weight. The second example (Fig. 7(b)) represents an initial arbitrary bone microstructure (trabecular) distributions onto which angular mechanical load are applied to observe the microstructure evolution for bone remodeling. Finally, case three (Fig. 7(c)) represents the same initial arbitrary bone microstructure (trabecular) distributions onto which a compressive triangular load is applied. Cases two and three are displacement fixed, only bone remodeling is allowed. Results are presented for both the macroscopic bone density evolution as well as mesoscopic trabeculae reorientation.
The presented theoretical numerical model was used on both the continuous and heterogeneous distributions. Results are presented on Fig. 8(a) and (b) for the continuous macroscopic geometry. They show a good correlation with 2D data available in the literature (see Carter et al., 1989 [43]). However, the development of a 3D model is able to provide more information since most full scale numerical models for bone reconstruction show mainly 2D numerical analyses. For any numerical models, identification of the parameters will be different between 2D and 3D analyses. On Fig. 8(b), we can observe clearly that bone density is not constant through the thickness which in turn will influence directly the bone density evolution. In addition, the model parameters being biological (cell activation), it is possible to extract localized information in order to propose a more specific multiscale, multiphysics approach for more patient specific analyses.

Macroscopic evolution from applied external mechanical boundary conditions.
Accounting for the local cell distribution in the macroscopic volume, we used the current model on an arbitrary oriented bone microstructure with two different loading conditions as defined on Fig. 7(b) and (c) to look more specifically at the bone microstructure evolutions as a function of the cells activation process. Results are presented on Figs 9 and 10 for the first and the second load conditions respectively. Bone density (in %), bone density gradient (in %/unit time) and mechanobiological stimulus are represented for different times, namely initial (or starting condition), intermediate (or during the evolution) and final equilibrium condition.
On the load case 1 (Fig. 9), no structural displacements are allowed, only bone remodeling. Bone density evolution (in %) shows that trabecular reorientation is mainly directed in the main principal stress axes (following main load directions) and that trabeculae thicknesses depend directly on the mechanical equilibrium. Local cells activation is observed through the variation of bone gradient on line two. It can be extracted and needs to be quantified for real medical application. The link between the second and third line shows more specifically the cell activation through the mechanobiological stimulus. The bone density variation is directly dependent on the localization of the calculated stimulus.

Mesoscopic (trabecular size) evolutions for the 45° load condition.
On Fig. 10, a similar analysis was run but using different external boundary conditions. A triangular compressive load was applied on the same initial bone microstructure. Similar observations can be made compared with the first load case. The trabeculae are mainly oriented in the principal load directions with decreasing size with lower load intensities. Localization of the cell activation depends on the calculated mechanobiological stimulus and needs to be quantified in the same manner.

Mesoscopic (trabecular size) evolution for the triangular load condition.
A comparison of the obtained results for triangular load is presented on Fig. 11 with data obtained on a phenomenological model from the literature [44].

Comparison of bone density distribution between current mechanobiological model and phenomenological model after applying triangular compressive mechanical load.
Here, good correlation is observed in terms of bone density distribution. Only, the current model enables a better quantification of the biological activations and a more patient dependent approach. The obtained bone microstructure is completely dependent on the mechanical biological parameters and can be quantified. Mechanical equilibrium is reached through biological cells activation and localization, and can be correlated with a physical reality.
Similar works, linking local mechanics with biology, were undertaken by Ruimermann [45] to extract specific patient bone remodeling processes for 3D micro-architectures but no links were made with real medical applications. The development of the current model will help in bridging these scales to be able to obtain a predictive macroscopic mechanobiological model for bone evolutions as a function of the local biological constituents.
We presented here some works related to the interests of using theoretical multiscale numerical models for the prediction of soft and hard biological tissue behaviors. These models showed the possibilities of extracting specifically the effects of the local constituents on macroscopic behaviors, be patient dependent and enable good predictions for long term evolutions. Theory, numerical modeling and experiments are well correlated and show the need of pursuing these studies. Further works are still needed to improve the models accuracy and more experimental data will help to provide better patient dependent predictions.
Conflict of interest
None to report.
