A multiscale data-driven approach for bone tissue biomechanics

A multiscale data-driven approach for bone tissue biomechanics J. Mora-Mac´ıasa, J. Ayensa-Jim´enezb, E. Reina-Romoc, M.H. Doweidarb, J. Dom´ınguezc, M. Doblar´eb, J.A. Sanz-Herrerac,∗ aUniversity of Huelva, Huelva, Spain bMechanical Engineering Department, University of Zaragoza, Spain Aragon Institute of Engineering Research (I3A), University of Zaragoza, Spain Centro de Investigaci´on Biom´edica en Red en Bioingenier´ıa, Biomateriales y Nanomedicina (CIBER-BBN), Spain cEscuela T´ecnica Superior de Ingenier´ıa, Universidad de Sevilla Abstract The data-driven methodology with application to continuum mechanics relies upon two main pillars: (i) experimental characterization of stress-strain pairs associated to different loading states, and (ii) numerical elaboration of the elasticity equations as an optimization (searching) algorithm using compatibility and equilibrium as constraints. The purpose of this work is to implement a multiscale data-driven approach using experimental data of cortical bone tissue at different scales. First, horse cortical bone samples are biaxially loaded and the strain fields are recorded over a region of interest using a digital image correlation technique. As a result, both microscopic strain fields and macroscopic strain states are obtained by a homogenization procedure, associated to macroscopic stress loading states which are considered uniform along the sample. This experimental outcome is here referred as ∗Corresponding author. Camino de los descubrimientos s/n, 41092 Seville, Spain. Tel.: +34 954 486079; fax: +34 954 487295 Email address: [email protected] (J.A. Sanz-Herrera) Preprint submitted to Elsevier June 22, 2020 This is an Accepted Manuscript of an article published by Elsevier in Computer Methods in Applied Mechanics and Engineering,Volume 368, 15 August 2020, available at: https://doi.org/10.1016/j.cma.2020.113136 //Copyright 2020 Elsevier. En idUS: Licencia Creative Commons CC BY-NC-ND a multiscale dataset. Second, the proposed multiscale data-driven methodology is implemented and analyzed in an example of application. Results are presented both in the macroscopic and microscopic scales, with the latter considered just as a post-process step in the formulation. The macroscopic results show non-smooth strain and stress patterns as a consequence of the tissue heterogeneity which suggest that a preassumed linear homogeneous orthotropic model may be inaccurate for bone tissue. Microscopic results show fluctuating strain fields –as a consequence of the interaction and evolution of the microconstituents– an order of magnitude higher than the averaged macroscopic solution, which evidences the need of a multiscale approach for the mechanical analysis of cortical bone, since the driving force of many biological bone processes is local at the microstructural level. Finally, the proposed multiscale data-driven technique may also be an adequate strategy for the solution of intractable large size multiscale FE2computational approaches since the solution at the microscale is obtained as a postprocessing. As a main conclusion, the proposed multiscale data-driven methodology is a useful alternative to overcome limitations in the continuum mechanical study of the bone tissue. This methodology may also be considered as a useful strategy for the analysis of additional biological or technological hierarchical multiscale materials. Keywords: Data-driven approach, Computational biomechanics, Experimental bone tissue mechanics, Numerical simulation, Multiscale analysis 2 1. Introduction 1.1. Data-Driven approach in computational mechanics and multiscale analysis Data science (DS) has disrupted in the last years in industry along a variety of diverse fields, such as marketing and e-commerce [1], social sciences [2], healthcare [3] or the internet of things [4]. DS has also impacted in continuum mechanics with classical model-based material constitutive laws being progressively replaced by measured discrete stress–strain couples, which provide a more accurate and realistic representation of the physical relationships. The data-driven (DD) formulation in continuum mechanics was presented a few years ago in the context of elastostatics [5] and extended to elastodynamics in [6], as an alternative to classical model-based approaches. Related problems of the DD methodology, such as the unavoidable noise coming from experimental testing campaigns, were addressed in [7] and in [8]. Moreover, [9] presented a theoretical description to fill and expand (upsampling) missing, incomplete and unreliable (noisy) datasets. This latter framework is used in this work in a particular case of interest. The DD mechanics scheme proceeds by searching the closest strain-stress couple in the measuring set, through the minimization of a certain metric –which may eventually include a stochastic definition [8]–. The optimization problem is then constrained to fulfill the equilibrium and compatibility equations, which are introduced by means of Lagrange multipliers [5]. An alternative formulation was presented in [10] introducing the concept of constitutive manifold, which embeds the stress-strain couples, while the found solution satisfies all the required equations. The description of the construc- 3 tion of such manifold is detailed in [11]. On the other hand, multiscale simulations in classical continuum mechanics, establishes the analysis at different spatial and temporal observation scales such that the response of the upper (coarse) scale is a direct consequence of the evolution of the lower (fine) scale. The purpose of multiscale analyses is two-fold: First, to establish less phenomenological models based on more fundamental principles and spatio-temporal resolutions. Second, to have access to the distribution of the variables along both the fine and coarse scales. This kind of multiscale approaches are especially well-suited to the mechanical analysis of heterogeneous materials, which contain an underlying microstructure such as natural (biological) or artificial (technological) materials. Different multiscale techniques have been implemented to account for plasticity and the subsequent microstructure evolution [12, 13], poroelasticity and microstructural fluid circulation [14] or micromechanical analysis of composite adhesives [15, 16] to cite a few. The multiscale approach has also been applied in the framework of bone tissue mechanics. In these applications, loads are transferred at the macroscopic organ scale whereas the driving stimuli of the internal biological evolution processes are local and microstructural [17, 18, 19]. Several examples for bone tissue engineering can be seen in [20, 21, 22] among other. The reader is addressed to [23] for a review of multiscale techniques applied to continuum models. Multiscale analysis based on finite elements are usually called FE2[24, 25] due to both the complexity and CPU time scale to the power of 2: the corresponding simulations have to be performed both at the macro and micro levels, this latter defined at each Gauss point of the macro finite element 4 mesh. Even though some remedies have been proposed in the literature to alleviate the CPU time demand [26], this kind of analyses is usually intractable with limited applicability to the practical use. The present work exploits the experimental data acquisition of the mechanical characterization of a bone specimen, available both at the tissue (macro) and microstructural (micro) levels, to propose a multiscale approach in the context of DD. It is clear that the feasibility of this DD technique relies on the availability of test data at both scales. However, once a proper macro and micro tissue characterization is performed, the solution at the microstructure is obtained as a postprocessing in the DD methodology, avoiding prohibitive CPU times. The experimental protocol used in this work, in the context of bone tissue mechanics, is reviewed next. 1.2. Experimental characterization of bone tissues DD technique is inherently linked to the experimental characterization of the material in hands, such as the bone tissue in this study, prior to DD simulations. An example of bone tissue is cortical bone, which is found in the diaphysis of long bones as an outer layer of the trabecular bone. At the macroscopic level, it looks like a compact and continuous solid. However, at the microscopic level, the cortical bone structure is very complex as seen in Fig. 1. Therefore, the mechanical response and properties measured at the macroscopic level of the cortical bone are influenced by its microarchitecture. This relation between micro and macroscale has been studied for decades. For example, [27] related the mechanical properties with the collagen fiber orientation. More recently, [28] showed that the nanoscale heterogeneity of the mechanical properties of the bone promotes energy dissipation. A 5 multiscale DD approach could take advantage from the data provided by experimental techniques which allows characterizing the cortical bone tissue at the microscopic level. Among the techniques used to characterize the cortical bone tissue at that microscopic level, histologies, and microscopy techniques allow characterizing biological aspects of the bone microarchitecture [27]. More recent techniques, such as microcomputed tomography, provide bone microarchitecture, and other parameters that may be related with mechanical properties such as the mineral density [29, 30]. Nanoindentation can also provide mechanical properties of a bone surface, which in combination with microscopic techniques or microcomputed tomography may define the bone heterogeneity and the spatial variations of its mechanical properties [28, 31]. However, none of these methods is able to measure the micromechanical response (stress and/or strain at the microscopic level) of the bone tissue during mechanical experiments. Digital image correlation (DIC) technique has shown advantages over the above mentioned methods to measure the microscopic strain field, of the order of tens of microns, in real time, in bone samples under different loading conditions. Consequently, the macroscopic strain field, of the order of millimeters, can be obtained by averaging along the domain of the region of interest, or the microstructural representative volume element (RVE). This technique consists in comparing the positions in a speckle pattern of reference, usually spray-painted at the loaded sample surface, where the images are taken. Refs. [32] and [33] used DIC to compute the strain field in a region of the surface of a mouse tibia, with a resolution of 125 µm distance among 6 the field strain points. Thompson et al. [34] also used DIC for determining the local mechanical conditions within the early bone callus providing strain fields of the same order of magnitude. Currently, micro-strain fields provided by DIC measurements have been integrated into experimental studies which combine different length scale measurements (in situ loading, digital image correlation, and synchrotron X-ray) in order to understand the deformation mechanisms within the bone tissue [35]. As a conclusion, DIC is a contrasted and validated technique to measure microscopic strain fields above the micron scale, and associated strain stress macroscopic field in bone tissue. Figure 1: Microarchitecture of cortical (compact) bone tissue [36]. 7 1.3. Objective and organization of the paper The main objective of the paper is to present a multiscale DD study in the framework of bone tissue mechanical behavior. According to the DD methodology, it includes an ad hoc mechanical characterization of the cortical bone tissue –through macroscopic stress–strain pairs and microscopic strain fields– by means of DIC measurements under different loads and boundary conditions. This information is used as a multiscale dataset (including a data completion technique) along a specific algorithmic development for the equations. The paper is then organized according to this scheme: Section 2 describes the bone tissue characterization using the DIC technique at the microstructural scale. On the other hand, the numerical treatment of the elasticity equations as well as the strategy to deal with the experimental results both at the macro- and micro-scales along the formulation, is also shown in section 2. Section 3 presents both tissue characterization results, as well as multiscale DD results in a selected example of application. The discussion is addressed in section 4, whereas some final conclusions are drawn at the end of the paper. 2. Materials and methods 2.1. Experimental setup, bone tissue characterization A sample of cortical bone was taken from the mid diaphysis of the femur of an adult horse. It was extracted in the longitudinal direction and then it was sawed and polished with carbide papers (P600 to P4000) and diamond slurry. The final dimension of the sample was 50 ×20 mm with a thickness of 4 mm (fig. 2a). The sample was subsequently embedded in epoxy resin to 8 fix it to the mechanical testing machine (fig. 2b). Then, it was spray-painted to generate a speckle pattern for the strain measurements as explained next. (a) (b) Figure 2: Sample preparation. (a) Sample of an adult horse femur after sawing and polishing; (b) the sample was subsequently embedded in epoxy resin and fixed to the threaded ends. The loading of the specimen was performed using a push-pull testing machine that allows applying macroscopic biaxially longitudinal (LF) as well as transversal (TF) compressive forces (see fig. 3). Both longitudinal and 9 to expand data coverage and to fill missing data or gaps in the sample. 2. It can be used appropriately as a field reconstructor. Indeed, a synthetic field may be computed from an existing dataset, but taking into account reliability and physical admissibility of the data. Here, the data-completion technique is used in these two ways. On the one hand, at the macroscale (M) level, pairs of strain-stress couples εMand σM are sampled from an obtained experimental setup as explained in the previous section. This give us an unstructured dataset E={(εM,i,σM,i)}i=1,...,N , where Nis the sample size. However, this amount of data Nmay not be large enough for some purposes. For instance, if there is no control in one of the sampled variables, such as in longitudinal strains and stresses in our experimental setup, data coverage may be poor. Besides, for computational purposes, the resolution of data is often not high enough, needing some tools to work, sometimes at the continuum level. The methodology presented in [9] is therefore used in order to obtain, from a unstructured dataset E, a new structured data-set satisfying ui=umin +i n(umax −umin), i= 1, . . . , n, nbeing the number of upsampled points. umin and umax are the minimal and maximal values reached by the variable u.u=εM,σM. The filling procedure is performed by the following strategy: for each stress state σM,∗, a new strain state εM,∗is computed by solving the minimization problem, εM,∗= argmin εM N X i=1 ϕdσM,i,σM,∗d2εM,i,εM(2) where dis the euclidean distance and ϕ:R+→[0; 1] is a suitable activation function. This problem is easily solved by formulating a linear system, as explained in [9]. 16 On the other hand, the data-completion technique was also used for the computation of the microscale (m) strain field. In this case, the aim is to obtain a microscale field, built from the data-set of possible microscale fields, but compatible with the macroscale pairs of strains and stresses. The approach here is similar, except for the fact that we must impose the constraint relating the macroscopic strain εMand the microscale strain field εm(y) as expressed in equation (1). As only the microscopic strains fields compose the multiscale dataset, these are the only data that may be used in the minimization procedure to obtain a microscale field. Therefore, the problem is here reduced to compute εm,∗(y) by solving the following minimization problem: εm,∗(y) = argmin εm(y) N X i=1 ϕd1 VZV εm,i(y)dV, εMd2εm(y),εm,i(y) subject to 1 VZV εm(y)dV =εM (3) Note that 1 VRVεm,i(y)dV is the macroscale strain associated to the microscale field εm,i(y), εMis the macroscopic strain whose microscale strain field is intended to be obtained and E={εm,i(y)}i=1,...,N. is the learning dataset. As the constraint is linear, this problem is easily solved following the two-steps strategy described in [9], where first a linear system is solved and then the solution is projected into a manifold, obtaining the final result. Of course, the accuracy of the result depends on the coverage of the dataset. For the data presented in this work, a Radial Basis Function (RBF), ϕ(u) = exp(−u2 ζ), was selected in Eqs. (2) and (3) as activation function, with ζthe spread parameter. 17 Finite element Gauss point Stress Strain Macrodata Microdata Figure 6: Numerical scheme of the multiscale DD implementation in a finite element framework. 18 2.2.2. Multiscale data-driven computational approach In this section, the multiscale DD computational approach and its numerical implementation are exposed in the framework of a finite element discretization. The implementation is an extension of the algorithm presented by [5], developed at small strains, and proceeds by searching (macroscopically) the closest stress-strain pair, i.e. σM,εM, to the multiscale dataset D≡σM,εM×εm(y)D. Note that superindices Mand mrefer to the macro and micro datasets, see fig. 6. That is, the dataset contains both the macroscopic stress-strain pairs, as well as the associated microscopic distribution (fluctuation) of the strain field εm D(y) along the microscopic RVE domain y∈ V (fig. 6). In fact, the macroscopic strains are obtained from the microscopic strain distributions in the RVE as stated in equation (1). The problem is established macroscopically, leading to an optimization problem that must be constrained to the equilibrium and compatibility (macroscopic) equations. This may be written as follows: min(σM,εM,σM D∈D,εM D(εm D)∈D)1 sd2(σM,σM D) + 1 d2(εM,εM D(εm D)) s.t. ∇ · σM=0 εM=1 2∇uM+∇TuM σM·nM=tMin Γt uM=¯ uMin Γu (4) with uMthe macroscopic displacement field, sand representative values of the stress and strain ranges in the test data, respectively. tMare the prescribed values of the macroscopic traction field at the Neumann boundary Γt, 19 associated to the normal vector nM, while ¯ uMare the prescribed macroscopic values of the displacement field at the Dirichlet boundary Γu.dis the metric defining the distance of the pair at hand from the dataset. This metric was defined as the Euclidean distance, d(x,xD) = sX i (xi−xi,D)2(5) Substitution of Eq. (5) into (4) yields, min(σM,εM,σM D∈D,εM D(εm D)∈D)n1 s σM−σM D  2+1  εM−εM D(εm D)  2o s.t. ∇ · σM=0 εM=1 2∇uM+∇TuM σM·nM=tMin Γt uM=¯ uMin Γu (6) Using the finite element (FE) discretization procedure over a single element, Eq. (6) yields, minσ(e),M ,ε(e),M ,σ(e),M D∈D,ε(e),M D(εm D)∈D1 s  σ(e),M −σ(e),M D   2+1   ε(e),M −ε(e),M D(εm D)   2 s.t. ε(e),M =B(e)·u(e),M Pkwk·B(e),T k·σ(e),M k=F(e),M (7) where B(e)is the shape functions gradient matrix and F(e),M the vector of macroscopic forces which accounts for the body and boundary forces. Subscript (e) denotes the discretized variables defined over an element. On the 20 other hand, wkare the corresponding weights of the integration (Gauss) points defined over the element. The first compatibility (discretized) constraint in (7) can be directly substituted back in the objective function. On the other hand, the second FE discretized constraint can be included via a Lagrange multiplier, finally resulting in the following expression: minσ(e),M ,ε(e),M ,σ(e),M D∈D,ε(e),M D(εm D)∈D1 s  σ(e),M −σ(e),M D   2+1   B(e)·u(e),M −ε(e),M D(εm D)   2+ +Pkwk·B(e),T k·σ(e),M k−F(e),M ·η(e),M  (8) with η(e),M being the discretized Lagrange multiplier defined over the element (e) of the macrostructure. Taking variations in (8), we get: δu(e),M =0→1 B(e)·u(e),M −ε(e),M D(εm D)=0(9a) δσ(e),M =0→1 sσ(e),M −σ(e),M D+B(e)·η(e),M =0(9b) δη(e),M =0→X k wk·B(e),T k·σ(e),M k−F(e),M =0(9c) After some algebraic manipulation of (9b) and substitution of (9c) in (9b) we obtain, Pkwk·B(e),T k·B(e) k·u(e),M =Pkwk·B(e),T k·ε(e),M D,k (εm D) Pkwk·B(e),T k·s·B(e) k·η(e),M =Pkwk·B(e),T k·σ(e),M D,k −F(e),M (10) Eq. (10) is then integrated over all finite elements of the mesh. Hence, ANel e=1 Pkwk·B(e),T k·B(e) k·uM=ANel e=1 Pkwk·B(e),T k·ε(e),M D,k (εm D) ANel e=1 Pkwk·B(e),T k·s·B(e) k·ηM=ANel e=1 Pkwk·B(e),T k·σ(e),M D,k −FM(11) 21 where Ais the FE assembly operator and Nel the number of elements of the FE mesh. Note that FMrefers to the global macroscopic force vector after assembly. The algorithmic implementation of the set of Eqs. (11) follows the scheme presented by [5], but here for multidimensional stress-strain pairs. The numerical implementation of the algorithm was performed in Matlab. The selection of the stress-strain pairs in the dataset is divided into normal and shear components. Assuming orthotropic symmetries in the bone material along the longitudinal and transversal directions, normal components of the pairs (σM,εM) are selected from the dataset (with independence on the values of shear compornents) as the closest distance to (σM,εM)D(following the Eulerian norm in (5)) using the Matlab function knnsearch. The corresponding shear stress and strain components are computed indirectly through principal components as follows: εM αβ =Qαi ·Qβj ·εM ij (12a) σM αβ =Qαi ·Qβj ·σM ij (12b) with Qαi the principal directions of the strain tensor εij. Indices i, j are referred to L, T (longitudinal and transversal) directions, and α, β to I, II principal directions. Then, the closest point of the dataset (σM I, σM II ;εM I, εM II )D to (σM I, σM II ;εM I, εM II ) is obtained using the knnsearch function. Finally, the stress and strain shear components are recovered as follows: εM LT =QLα ·QT β ·(εM αβ)D(13a) σM LT =QLα ·QT β ·(σM αβ)D(13b) 22 Note that obtaining the stress and strain shear components as explained above is equivalent to assume isotropy to obtain the shear compoenents. We remark that we use Eqs. (13) only to compute the shear components of the stress and strain tensors but not the normal components. 3. Results 3.1. Experimental results Fig. 7 shows the mean values of the longitudinal and transversal strains obtained for different applied loads. It may be observed that the minimum longitudinal strain is reached for the maximum longitudinal compression (LF =-1500 N) and the minimum transversal compression (TF = 0 N), since the latter increases the longitudinal strain due to Poisson effect. An opposite trend may be observed in the case of the transversal strain, whose minimum values are obtained for the minimum longitudinal compression (LF = 0 N) and the maximum transversal compression (TF = 600 N). Using a standard orthotropic finite element model of the experimental setup, it was corroborated that macroscopic stress and strain fields were uniform across the ROI. Indeed, we concluded that the stress can be estimated in the ROI as the applied force over the compressive surface. Assuming linear elasticity and looking at the longitudinal strain and force, it is possible to estimate an elastic modulus of the tested cortical bone sample by dividing the maximum applied longitudinal stress, 18.75 MPa (compression area 4×20 mm) by the maximum longitudinal strain, around 0.1%. The obtained estimation of 18.75 GPa is within the order of magnitude of the values published in previous studies for cortical bone [31]. 23 Figure 7: Macroscopic load states applied over the sample and mean obtained values for longitudinal (εM L) and transversal (εM T) strains. Each of the macroscopic load states represented in fig. 7 provides a dataset of 192 ×176 points, with the associated microscopic strain field. It is used as an input multiscale dataset of the computational analysis. As an example, the microscopic strain fields for the experimental macroscopic load states, as detailed in table 1, are represented in figs. 8 and 9. 3.2. Completion of the dataset The results regarding the measured dataset, both macroscopic and microscopic data, are expanded by means of the mathematical approach presented in section 2.2.1 (see [9]). The macroscopic upsampled data are presented in fig. 10 using a sampling size of n= 100 ×100 and a spread parameter (force) ζ= 80 N. The upsampled macroscopic dataset is presented over the 24 (1) (2) (3) (4) Figure 8: Experimentally obtained microscopic longitudinal strain fields (εm L) for the experimental macroscopic load states detailed in table 1 (1, 2, 3 and 4). 25 Figure 14: Estimated (reconstructed) microscopic strain fields, by means of the data completion technique, versus real ones measured by means of DIC. Randomly selected macroscopic strain state (εM T,D, εM L,D) indicated in the figure. 32 MPa, respectively, and to prescribed normal displacements as described in fig. 15. The FE mesh was composed of 1240 triangular linear elements. The used upsampled dataset micro-macro dataset includes n= 500 ×500 longitudinal and transversal points. The values sand in Eq. (4) are selected as the experimental range of stress and strain in the dataset, respectively. We obtained (data not shown) that different selected values (although keeping the same order of magnitude) converged to the same results. The macroscopic longitudinal and transversal components of the displacement field are shown in figure 16. Moreover, the macroscopic longitudinal and transversal components of the strain and stress tensors are shown in figs. 17 and 18, respectively. 3.4. Computational multiscale DD results Multiscale results are presented along the microstructure, i.e. RVE, at three selected points in the macroscale. The distribution of the microscopic strain field in the RVE represents the fluctuation of the strain field as a consequence of the heterogeneity of the domain at this scale. Figs. 19, 20 and 21 show three points in the macroscale with the associated macroscopic strain fields and their underlying microstructural distribution of the corresponding strain fields in the RVE. 4. Discussion We may set the beginning of the application of DD techniques to continuum mechanics just a few years ago [5]. Since then, several papers have followed this research direction. Even though DD needs experimental datasets as input, most of the works use synthetically generated artificial data. This 33 LF TF 10 cm T 10 cm 5 cm 2.5 cm 45o L Figure 15: Example of application for multiscale DD simulations. 34 Figure 16: Macroscopic displacement field [m] of the example of application shown in fig. 15. Left: macroscopic longitudinal displacement. Right: macroscopic transversal displacement. 35 Figure 17: Macroscopic strain field [–] of the example of application shown in fig. 15. Left: macroscopic longitudinal strain. Right: macroscopic transversal strain. Figure 18: Macroscopic stress field [MPa] of the example of application shown in fig. 15. Left: macroscopic longitudinal stress. Right: macroscopic transversal stress. 36 Macro Micro Figure 19: Multiscale DD results. Top: Macroscopic strain field [–] of the example of application shown in fig. 15. Bottom: Microscopic strain field in the RVE for the highlighted point of the macroscale. Left: longitudinal strain component. Right: transversal strain component. 37 Macro Micro Figure 20: Multiscale DD results. Top: Macroscopic strain field [–] of the example of application shown in fig. 15. Bottom: Microscopic strain field in the RVE for the highlighted point of the macroscale. Left: longitudinal strain component. Right: transversal strain component. 38 Macro Micro Figure 21: Multiscale DD results. Top: Macroscopic strain field [–] of the example of application shown in fig. 15. Bottom: Microscopic strain field in the RVE for the highlighted point of the macroscale. Left: longitudinal strain component. Right: transversal strain component. 39 strategy neglects important problems and limitations that the DD methodology has to face when dealing with experimental setups in realistic applications. We are referring to these limitations below in this section. The present work uses a real experimental set of data at multiple scales to characterize a piece of cortical bone tissue, a material with a high interest in the biomedical community. The mechanical characterization was performed by means of the DIC technique similar to other works in this context [32, 33, 35]. The conducted experimental setup allowed obtaining information both at the micro- and macro-scales in a multiscale fashion. According to figs. 8 and 9, the fluctuation of the microscopic strain fields as a consequence of the heterogeneity at this scale –due to the porous internal architectures with heterogeneous mineral content (see fig. 5a)– justifies a multiscale approach of the cortical bone tissue. The multiscale dataset was then expanded (upsampled) following a previous approach proposed by the authors [9]. This upsampling procedure is often necessary in real situations since the resolution of the experimental stress or strain driven setups may be coarser than the one required to get an adequate resolution in DD simulations [42]. The proposed multiscale DD approach was tested in a 2D plane stress geometry. It must be clearly remarked that DD simulations are inherently linked to the availability of data. Since the experimental characterization was established in a 2D setup (biaxially loaded states), simulations are therefore restricted to 2D cases. Note that, in general terms, 3D mechanical characterization is not straightforward. This fact can be identified as a limitation in the DD methodology, directly derived from our current experimental ca- 40 pability. The macroscopic strain and stress fields, shown in figs. 17 and 18, respectively, provide an unexpected non-smooth solution versus a smooth one which could be achieved using a classical model-based linear orthotropic approach for the constitutive response of the cortical bone tissue [43, 44, 45]. This observed behavior is a direct consequence of the measured non-smooth pattern in the dataset (Fig. 10). Therefore, we hypothesize that this nonsmooth characteristic is a result of the microstructural heterogeneity of the material in hands. This hypothesis is enforced by the presence of the internal microstructure and fluctuation of the field variables, as presented in Figs. 19, 20 and 21. This kind of macroscopic non-smooth pattern for strain and stress fields is typical of multiscale simulations which consider the internal evolution of the microstructure as can be seen in Refs. [23, 46, 47], amongst others. Moreover, these non-smooth patterns suggest that DD approaches are useful and give value at certain preassumed scenarios, or a priori wellknown behaviors such as the linear orthotropic model-based approach of the bone tissue. Even though the substitution of well-accepted constitutive laws by data may be considered a waste of knowledge [48], we have demonstrated in this work that this is not the case for the cortical bone tissue. Multiscale results give access to the microscopic strain distribution in the RVE. Figs. 19, 20 and 21 show that the peak fluctuations of the strain fields in the finer (microscopic) scale is an order of magnitude higher than the averaged (homogenized) strain fields. Since many of the physiological processes in the bone tissue are microstructural and local, including microcracking, bone fracture, and bone microarchitectural disorders (as menopause disease); information of this mechanical variable at the microscopic level is 41 [32] A. Carriero, L. Abela, A.A. Pitsillides, S.J. Shefelbine, Ex vivo determination of bone tissue strains for an in vivo mouse tibial loading model, J. Biomech. 47(10) (2014) 2490–2497. [33] P. Sztefek, M. Vanleene, R. Olsson, R. Collinson, A.A. Pitsillides, S. Shefelbine, Using digital image correlation to determine bone surface strains during loading and after adaptation of the mouse tibia, J. Biomech. 43(4) (2010) 599–605. [34] M.S. Thompson, H. Schell, J. Lienau, G.N. Duda, Digital image correlation: a technique for determining local mechanical conditions within early bone callus, Med. Eng. Phys. 29(7) (2007) 820–823. [35] A. Gustafsson, N. Mathavan, M.J. Turunen, J. Engqvist, H. Khayyeri, S.A. Hall, H. Isaksson, Linking multiscale deformation to microstructure in cortical bone using in situ loading, digital image correlation and synchrotron X-ray scattering, Acta Biomater. 69 (2018) 323–331. [36] E.N. Marieb, Bone skeletal tissues Part A. Human anatomy and physiology. Sixth edition. Pearson education Inc. (2004). [37] J.D. Currey. Bones: Structure and Mechanics. Princeton University Press Editions (2006). [38] J. V´azquez, Efecto de las Tensiones Residuales en la Fatiga por Fretting, PhD Thesis University of Seville, 2010. [39] M.A. Sutton, S.R. McNeill, J.D. Helm, Y.J. Chao, Advances in twodimensional and three-dimensional computer vision, In: Rastogi PK, ed- 48 itor. Photomechanics, Topics Appl. Phys. Springer-Verlag Berlin Heidelberg, 77, 323–372, 2000. [40] M.A. Sutton, J.J. Orteu, H. Schreier, Digital Image Correlation (Chapter 5), In: Sutton M.A., Orteu J.J. and Schreier H., Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts Theory and Applications, pp. 81–116, 2009. [41] P.M. Suquet, Elements of homogenization for inelastic solid mechanics. Trends and applications of pure mathematics to mechanics E. Sanchez- Palencia, A. Zaoui (Eds.), Homogenization techniques for composite media, Lecture notes in physics, Springer, Berlin, pp. 193–278, 1987. [42] R. Eggersmann, T. Kirchdoerfer, S. Reese, L. Stainier, M. Ortiz, Model- Free Data-Driven inelasticity, Comput. Methods Appl. Mech. Engrg. 350 (2019) 81-99. [43] D. Reilly, A. Burstein, The elastic and ultimate properties of compact bone tissue, J. Biomech. 8 (1975) 393–405. [44] W.R. Taylor, E. Roland, H. Ploeg, D. Hertig, R. Klabunde, M.D. Warner, M.C. Hobatho, L. Rakotomanana, S.E. Clift, Determination of orthotropic bone elastic constants using FEA and modal analysis, J. Biomech. 35 (2002) 767–773. [45] S. Bernard, Q. Grimal, P. Laugier, Accurate measurement of cortical bone elasticity tensor with resonant ultrasound spectroscopy, J. Mech. Behav. Biomed. Mater. 18 (2013) 12–19. 49 [46] H-W. Zhang, J-K. Wu, J. Lu, Z-G. Fu, Extended multiscale finite element method for mechanical analysis of heterogeneous materials. Acta Mech. Sin. 26 (2010) 899–920. [47] K. Matous, M.G.D. Geers, V.G. Kouznetsova, A. Gillman, A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials. J. Comput. Physics 330 (2017) 192–220. [48] D. Gonzalez, F. Chinesta, E. Cueto, Learning corrections for hyperelastic models from data, Front. Mat. 6 (2019) 14. [49] H. Isaksson, O. Comas, C.C. van Donkelaar, J. Mediavilla, W. Wilson, R. Huiskes, K. Ito, Bone regeneration during distraction osteogenesis: mechano-regulation by shear strain and fluid velocity, J. Biomech. 40 (2007) 2002–2011. [50] A. Leygue, M. Coret, J. R´ethor´e, L. Stainier, E. Verron, Data-based derivation of material response, Comput. Methods Appl. Mech. Engrg. 331 (2018) 184–196. 50 Appendix: Computational cost estimation of the presented multiscale DD methodology and comparison to standard FEM2algo- rithms In the presented DD multiscale study, the computation of the microscopic fields is carried out as a post-process as discussed in the paper. Therefore, the computational cost of our study is one of the DD techniques. Next, we discuss the computational cost of a DD technique versus a multiscale one in a generic framework. The DD algorithm turns into an iterative solver solution as exposed in Kirchdoerfer and Ortiz [5]. The cost of each iteration is one of a FEM solver for a single load state. Hence, the DD computational cost is the following: DD ∼$(ND)·O(N)·Niters (14) O(N) is the computational cost of a FEM iteration and Nis the size of the FEM macroscopic mesh. $(ND)·O(N) is the computational cost of a FEM iteration corrected by the searching algorithm of the dataset at each iterations. This cost is dependent on the size of the dataset ND. In our case, $has been estimated in 2 for a computer with core processor i7 @3.4 GHz and RAM 32 GB and using the Matlab searching algorithm knnsearch. Niters is the number of iterations which is of the order of 100-1000 in a DD approach [7]. On the other hand, the cost of a FE2(generic nonlinear) multiscale (MS) approach is [24], MS ∼NLiters ·nliters ·NG·O(n) (15) 51 where NGis the Gauss points of the macroscopic mesh. NLiters and nliters are the macroscopic and microscopic nonlinear iterations (Newton-Raphson for example). O(n) is the computational cost of a FEM iteration and nis the size of the FEM microscopic mesh. Therefore, the comnputational cost of the proposed DD multiscale method is competitive versus multiscale numerical schemes if NLiters ·nliters ·NG·O(n)>> $(ND)·O(N)·Niters (16) In terms of order of magnitude NG∼N,NLiters ∼10, nliters ∼10, Niters ∼ 2000 in a conservative case scenario. Then Eq. (16) yields, N·O(n)>> 20 ·O(N) (17) Moreover, assuming that the macro and micro meshes nand Nare comparable in a conservative case (in fact the microstructure mesh is usually more detailed in a multiscale approach), Eq. (17) yields: N >> 20 (18) Eq. (18) is usually satisfied. 52