Multiscale Mathematical Modeling in Dental Tissue Engineering: Toward Computer-Aided Design of a Regenerative System Based on Hydroxyapatite Granules, Focussing on Early and Mid-Term Stiffness Recovery

We here explore for the very first time how an advanced multiscale mathematical modeling approach may support the design of a provenly successful tissue engineering concept for mandibular bone. The latter employs double-porous, potentially cracked, single millimeter-sized granules packed into an overall conglomerate-type scaffold material, which is then gradually penetrated and partially replaced by newly grown bone tissue. During this process, the newly developing scaffold-bone compound needs to attain the stiffness of mandibular bone under normal physiological conditions. In this context, the question arises how the compound stiffness is driven by the key design parameters of the tissue engineering system: macroporosity, crack density, as well as scaffold resorption/bone formation rates. We here tackle this question by combining the latest state-of-the-art mathematical modeling techniques in the field of multiscale micromechanics, into an unprecedented suite of highly efficient, semi-analytically defined computation steps resolving several levels of hierarchical organization, from the millimeter- down to the nanometer-scale. This includes several types of homogenization schemes, namely such for porous polycrystals with elongated solid elements, for cracked matrix-inclusion composites, as well as for assemblies of coated spherical compounds. Together with the experimentally known stiffnesses of hydroxyapatite crystals and mandibular bone tissue, the new mathematical model suggests that early stiffness recovery (i.e., within several weeks) requires total avoidance of microcracks in the hydroxyapatite scaffolds, while mid-term stiffness recovery (i.e., within several months) is additionally promoted by provision of small granule sizes, in combination with high bone formation and low scaffold resorption rates.


INTRODUCTION
The importance of mathematical modeling in dentistry and related fields has steadily increased over the last decades. Thereby, the most popular examples concern Finite Element models of the mandibular system, dating back to at least the early 1990s (Hart et al., 1992;Koritoth and Versluis, 1997). While these models have been continuously improved over recent years, the proper choice of mechanical material properties as key model input parameters has gained particular attention. It has become more and more customary to derive these parameters directly from computed tomography (CT) images of the investigated organs, be it through a more heuristic, regression analysis-based approach (van Ruijven et al., 2007), or based on the combination of X-ray physics and continuum micromechanics concepts (Hellmich et al., 2008). The latter approach explicitly considers the hierarchical organization of bone down to the microscopic scales of cellular activity, i.e., to those which are in the very focus of modern bioengineering approaches, including regenerative medicine strategies. In the present contribution, we wish to extend the application range of dental mathematical modeling, from the mechanics of standard mandibular systems, to latest developments in modern bioengineering approaches. In more detail, we here explore for the very first time how an advanced multiscale mathematical modeling approach may support the design of a provenly successful tissue engineering concept for regenerating large bone defects in the human mandible (Komlev et al., 2002(Komlev et al., , 2003. The latter concept employs double-porous, potentially cracked, single millimeter-sized granules packed into an overall scaffold material, which is then gradually penetrated and partially replaced by newly grown bone tissue.
The granules themselves, exhibiting diameters from a few hundred micrometers to one or two millimeters, result from a processing route based on the effect of immiscible fluids (Komlev et al., 2002(Komlev et al., , 2003. Key morphological features of these granules are seen in the left-hand column of Figure 1: Firstly, they contain pores of two different characteristic lengths: small pores, with a characteristic length of one to a few micrometersthese pores are termed "micropores" hereafter; and large pores, with a characteristic length of several hundred micrometersthese pores are termed "mesopores" hereafter. A composite of randomly oriented hydroxyapatite crystals and the micropores constitutes the "base material" making up the granules. Upon increasing the observation scale by several orders of magnitude, one can discern not only the aforementioned mesopores, but also cracks which pervade the individual granules. Finally, the scaffold material is made up of the above described granules, with pore space in-between-due to the characteristic length of these pores, which is of the order of the granule diameters, these pores are termed "macropores" throughout the present paper. It is in these macropores, where the regeneration process starts, i.e., where new bone tissue is formed after implantation of the scaffold systems. During this regeneration process, the newly developing scaffold-bone compound needs to attain the stiffness of mandibular bone under normal physiological conditions. In this context, the question arises how the compound stiffness is driven by the key design parameters of the tissue engineering system: macroporosity, crack density, as well as scaffold resorption/bone formation rates. We here tackle this question by combining the latest state-of-the-art mathematical modeling techniques in the field of multiscale micromechanics, as reviewed in Section 2.1, into an unprecedented suite of highly efficient semi-analytically defined computation steps resolving several levels of hierarchical organization, from the millimeter-down to the nanometer-scale. This includes several types of homogenization schemes, namely such for porous polycrystals with elongated solid elements, for cracked matrixinclusion composites, as well as for assemblies of coated spherical compounds; described in great detail in Sections 2.2-2.5. These mathematical developments allow for first-ever insights into the mechanical functioning of the investigated tissue engineering system, including the role the aforementioned design parameters, as is documented in Section 3, before the paper finds its conclusion in Section 4.

Consideration of Material Hierarchy -Introduction of Representative Volume Elements (RVEs)
In recent years, continuum micromechanics-based homogenization theory (Hill, 1963;Suquet, 1997;Zaoui, 1997Zaoui, , 2002 turned out as particularly suitable means for the mathematical modeling of the mechanical behavior of complex hierarchical material systems found in biology and biomedical engineering (Hellmich and Ulm, 2002;Hellmich et al., 2004a;Fritsch and Hellmich, 2007;Bertrand and Hellmich, 2009;Fritsch et al., 2009;Hamed et al., 2010;Fritsch et al., 2013;Scheiner et al., 2016). In this context, a material is understood as a microheterogeneous body filling a macro-homogeneous representative volume element (RVE) with characteristic length ℓ, ℓ ≫ d, d standing for the characteristic length of inhomogeneities within the RVE, and ℓ ≪ L, L standing for the characteristic lengths of geometry or loading of a structure built up by the material defined on the RVE. Notably, for achieving results characterized by a quite good accuracy of about 5%, d and ℓ need to be separated by not more than a factor of 2 (Drugan and Willis, 1996), while ℓ and L need to be separated by a factor of 5 to 50 (Kohlhauser and Hellmich, 2013).
In general, the microstructure within such an RVE is so complicated that it cannot be described in complete detail. Therefore, the microstructural description within an RVE is restricted to the choice of quasi-homogeneous subdomains (called material phases), which are characterized by the following properties: (i) their shapes, (ii) their volume fractions within the RVE, (iii) their mechanical properties, and (iv) their mechanical interactions. Based on these characteristics, one can FIGURE 1 | Investigated material system. Three-level representation of the hydroxyapatite-based granular biomaterial (column on the right-hand side), following the morphological features found in images on different observation scales (column on the left-hand side); the depicted images have been acquired by means of scanning electron microscopy (hierarchical level I) and µCT imaging techniques (hierarchical levels II and III). then derive the homogenized (upscaled) behavior of the material on the observation scale of the RVE, i.e., the relation between homogeneous deformations acting on the boundary of the RVE and resulting macroscopic (average) stresses. If a single material phase is micro-heterogeneous itself, its mechanical behavior can be estimated by means of introducing RVEs within this phase, with characteristic lengths ℓ 1 ≤ d, comprising again inhomogeneities with characteristic length d 1 ≪ ℓ 1 , and so on. Such an approach is referred to as multi-step homogenization. At sufficiently low observation scales, it may provide "universal" phase properties, i.e., properties which are invariant throughout an entire material class, such as all bone tissues occurring in vertebrates (Fritsch and Hellmich, 2007). For the material system investigated herein, i.e., for assemblies of double-porous hydroxyapatite granules with bone tissue optionally coating the granules, the following suite of RVEs, with increasing characteristic lengths, is introduced: • On hierarchical level I, a microporous, overall isotropic, hydroxyapatite polycrystal emerges, see the bottom of Figure 1, showing a scanning electron micrograph of the polycrystalline microstructure on the left-hand side, and a two-dimensional schematical representation of the (actually three-dimensional) RVE I on the right-hand side: The latter is composed of spherical micropores (with volume fraction φ polyHA micro , also called polycrystalline microporosity). These micropores interact mutually with randomly oriented cylindrical hydroxyapatite crystals, with volume fraction f polyHA HA = 1 − φ polyHA micro . The microporosity typically amounts to φ polyHA micro = 0.445 (Dejaco et al., 2012). The characteristic length of the polycrystalline RVE I is in the order of 10 µm, see the bottom line of Figure 1.
• On hierarchical level II, a mesoporous, cracked matrix material makes up the individual granules, see the center of Figure 1, showing a micro-computed tomography (µCT) image of the microstructure within a granule on the left-hand side, and a two-dimensional schematical representation of the (actually three-dimensional) RVE II on the right-hand side: Namely, penny-shaped cracks, with vanishing volume fraction, while being quantified in number and size through the crack density parameter ǫ, see Equation (2), and spherical mesopores, with volume fraction φ gran meso , are embedded in the polycrystal matrix with properties arising from the structure of RVE I. Within RVE II, the latter matrix fills the volume fraction f gran polyHA = 1 − φ gran meso . The mesoporosity typically amounts to φ gran meso = 0.189 (Dejaco et al., 2012). The characteristic length of RVE II is in the order of 1 mm, see the middle row in Figure 1.
• On hierarchical level III, a macroporous conglomerate material consisting of mesoporous, cracked hydroxapatite granules and newly grown bone tissue emerges, see the top of Figure 1: granules with the stiffness of RVE II described above and filling volume fraction f congl gran , are surrounded by layers of newly grown bone tissue, with volume fraction f congl bone and stiffness derived from the ultrasonic tests of Ashman and van Buskirk (1987). These coated spherical elements are assembled, in mutual contact, to a granular conglomerate with macropores, with volume fraction φ congl macro , in-between.
At the time of granule implantation, no bone tissue has been formed yet, and this initial configuration is characterized by f congl bone = 0.

Stiffness Tensor Homogenization at Hierarchical Level I: Elasticity of Porous Hydroxyapatite Polycrystal
The elastic behavior of interpenetrating non-spherical crystals with pores in-between can be particularly well represented by the self-consistent stiffness estimate introduced by Fritsch et al. (2006), where infinitely many solid phases oriented in all space directions as well as one spherical pore phase are embedded in a matrix with zero volume fraction and with the stiffness of the homogenized material itself. The corresponding stiffness tensor of the water-filled porous polycrystal at hierarchical level I of Figure 1 reads as where φ polyHA micro and f polyHA HA are the volume fractions of the micropores and the hydroxyapatite needles; C HA and C microφ , respectively, are the fourth-order stiffness tensors of the hydroxyapatite crystals and of the micropores, respectively; ϑ and ϕ are the Euler angles quantifying the orientations of the hydroxyapatite crystals; P polyHA cyl (ϑ, ϕ) and P polyHA sph , respectively, are the fourth-order Hill (or morphology) tensors related to cylindrical and spherical inclusions, respectively, embedded in a matrix made up of the microporous hydroxyapatite polycrystal; and I is the fourth-order unit tensor, the components of which are defined via the Kronecker delta δ ij (δ ij = 1 if i = j and δ ij = 0 if i = 1), namely I ijkl = 1/2(δ ik δ jl + δ il δ jk ). The double integrals in Equation (1) can be evaluated based on Stroud's integration equations (Stroud, 1971;Pichler et al., 2009). Details regarding the computation of the Hill tensors P polyHA cyl (ϑ, ϕ) and P polyHA sph are presented in the Supplementary Material of this paper.
Numerical evaluation of Equation (1) requires knowledge of stiffness tensors C HA and C microφ . The isotropic stiffness of hydroxyapatite is known from the experiments performed by Katz and co-workers (Katz and Ukraincik, 1971;Gilmore and Katz, 1982), yielding a Young's modulus of E HA = 114 GPa, and a Poisson's ratio of ν HA = 0.27, see also (Hellmich and Ulm, 2002;Hellmich et al., 2004b). We here consider the case where the pore fluid is free to leave the microporosity upon loading of RVE I. This relates to so-called drained conditions, with C microφ = 0 (Thompson and Willis, 1991).

Stiffness Tensor Homogenization at Hierarchical Level II: Elasticity of Cracked Mesoporous Granule Material
Matrix-inclusion composites are preferentially represented by a so-called Mori-Tanaka-type morphology (Mori and Tanaka, 1973;Benveniste, 1987). At hierarchical level II, two types of inclusions are embedded into a matrix made of the porous polycrystal with a stiffness resulting from the homogenization scheme of Section 2.2: (i) spherical pores (the volume fraction of which is the mesoporosity φ gran meso ), and (ii) penny-shaped (non-frictional, open) cracks oriented in all space directions. The latter fill an only negligible volume fraction, so that their amount is quantified through their number per volume Frontiers in Physiology | www.frontiersin.org N and their radius r cr , combined into the so-called crack density parameter according to Budianksy and O'Connell (1976): Based on the works of Deudé et al. (2002) and Dormieux et al. (2004), the corresponding stiffness estimate of the mesoporous, pre-cracked granule material reads as where f gran polyHA and φ gran meso are the volume fractions of the matrix made up of the microporous hydroxyapatite polycrystal and of the mesopores; C polyHA is the fourth-order stiffness tensor of the microporous hydroxyapatite polycrystal matrix (see Section 2.2); C mesoφ is the fourth-order stiffness tensor of the mesopores, defined analogously to the stiffness tensor of the micropores (i.e., drained), and P polyHA sph is the Hill tensor for spherical inclusions embedded in the isotropic microporous hydroxyapatite polycrystal matrix, see the Supplementary Material for details. Tensor , also occurring in Equation (3), is defined via the Poisson's ratio of the microporous hydroxyapatite polycrystal, ν polyHA , see (Dormieux et al., 2004).

Microstress and Microstrain Fields in Bone Tissue-Coated Hydroxyapatite Granules and in the Macropores -Matrix-Inclusion Problems of the Hervé-Zaoui and of the Eshelby Type
Having reviewed and confirmed, respectively, the reliability of advanced self-consistent and Mori-Tanaka homogenization schemes at hierarchical level I (in Section 2.2) and at hierarchical level II (in Section 2.3), we now assign corresponding stiffness properties to individual granules, which we surround by increasingly thick coatings of newly formed bone tissue, and which we then pile up to a conglomerate of bone tissue-coated granules, see the granule and bone phases in RVE III shown in Figure 1. This goes beyond extending the classical self-consistent and Mori-Tanaka estimates, but requires the consideration of shell-like morphologies as pioneered by Hervé and Zaoui (1993); and this adaption motivates the following further developments: While the phase strains at hierarchical levels I and II were all estimated from homogenous inclusion strains of the underlying Eshelby-problem consisting of a spherical or cylindrical inclusion being embedded either into a real matrix, or into a fictitious matrix with the stiffness of the homogenized material, only the macropore phase exhibits this homogeneous property as concerns the RVE of hierarchical level III. In turn, as for the coated granule phase, Hervé-Zaoui's matrixcoated inclusion problem is considered (see Figure 2). This problem is characterized by the following boundary and field conditions: • Homogeneous strains at any location which is infinitely far from the inclusion center: where x is the position vector, ξ is the displacement field, and E 0 is the field of uniform displacements to which the auxiliary matrix with the properties of the overall conglomerate is subjected infinitely far from the inclusion.
Homogenization of the stiffness of the macro-porous granular scaffold material considering bone ingrowth, with spherical granules coated by bone tissue, and embedded in a polycrystal-type composite material consisting of the bone-coated granules and macropores.
Frontiers in Physiology | www.frontiersin.org ∂σ ϑr ∂r + 1 r • Kinematic relation: where ∇ s represents the symmetric gradient operator. Evaluation of Equation (10) in spherical coordinates yields (Salençon, 2001) The nature of the coated inclusion is reflected by stiffness properties being defined as functions of the radius measured from the center of that inclusion: with r gran and r bone denoting the radii of the individual granule and of the outer surface of the bone coating; with σ (r) and ε(r) denoting the stresses and strains prevailing in the Hervé-Zaoui coated inclusion problem; with C gran as the isotropic stiffness tensor of the granule material, determined according to Equation (3); with C bone as the stiffness tensor of newly grown bone tissue; and with C congl as the stiffness of the overall macroporous scaffold-bone compound of hierarchical level III. Due to the random orientation of the morphological features on all considered observation scales (hierarchical levels I -III), the overall conglomerate stiffness is isotropic, hence it is defined as with k congl and µ congl as the bulk and shear modulus of the scaffold-bone conglomerate on the macroscopic observation scale (hierarchical level III). As regards the stiffness tensor of newly formed bone tissue, C bone , we adapt the strategy described in Bertrand and Hellmich (2009) for the current purpose: We start with the ultrasoundbased tissue stiffness reported by Ashman and van Buskirk (1987), amounting to where the directions 1, 2, and 3 refer to the radial, circumferential, and axial directions of the orthotropic bone tissue material. The axial direction is aligned with the collagen fibril orientation, and during bone regeneration, the latter follows given morphological features occurring under normal physiological conditions. Since the spherical granular features deviate from the aforementioned physiological situation, it is very probable that the collagen fibrils are oriented along the tangent planes to the granule spheres, while not having any preferred orientation within these planes. As a first-order approximation of this situation, we let the stiffness tensor given by Equation (16) rotate first around axis 3 and average over all corresponding results (this leads to a transversely isotropic stiffness tensor with the isotropic plane coinciding with the 1-2 plane), and we then let the latter stiffness tensor rotate about an axis orthogonal to axis 3, which in turn leads to yet another transversely isotropic stiffness tensor with the plane of isotropy now including the axis 3. With respect to a spherical coordinate system attached to the granule, see Figure 2, the isotropic plane of the newly grown bone tissue coincides with the e ϑ -e ϕ -plane, so that C bone can be given as 15.90 9.00 9.00 0 0 0 9.00 21.74 10.70 0 0 0 9.00 10.70 21.74 0 0 0 0 0 0 11.04 0 0 GPa . (17) Given the isotropic nature of the overall scaffold-bone compound material, our interest in auxiliary homogenous strain fields E 0 prescribed at the infinite boundary of the Hervé-Zaoui matrix according to Equation (6) can be restricted to purely volumetric and purely deviatoric cases; as will become fully evident during the bulk and shear modulus homogenization described in Section 2.5. For the volumetric load case, the boundary condition defined by Equation (6) is restricted to where I is the second-order unit tensor, and e r is the base vector in radial direction, see Figure 2. Spherical symmetry of the remote loading conditions according to Equation (18) implies spherical symmetry of the resulting displacement field, with the corresponding spherically symmetric strain field reading as (20) Insertion of this strain field into the elastic laws given by Equations (12) -(14), and of the corresponding result into the equilibrium condition given by Equations (6) -(9) yields When introducing parameter n i , the general solution of the ordinary differential equation (21) reads as For the definitions of material parameters Ŵ k i,j (j = 1, 2), see the Supplementary Material. Evaluating Equation (22) for an isotropic phase, where C iso i,rrrr = C iso i,ϑϑϑϑ = k i + 4 3 µ i , and C iso i,ϑϑϕϕ = C iso i,rrϑϑ = k i − 2 3 µ i , one can easily see that n i = 3/2; this is the case for i = gran, scaff. For the bone tissue stiffness according the Equation (17), n bone amounts to 1.79. Inserting the general solution of the displacement field, Equation (23), into the kinematic relation, Equation (11), yields the corresponding strain field, given by its components in spherical coordinates, The stress field resulting from hydrostatic deformations follows from insertion of Equations (24) -(26) into the constitutive relations, Equations (12) - (14), as For the definitions of material constants M i , N i , O i , and P i , see the Supplementary Material. As regards the deviatoric loading, we prescribe a purely deviatoric (pure shear) deformation E d,0 at the infinitely remote boundary of the domain depicted in Figure 2, where e x and e y are the base vectors of a Cartesian coordinate system with its origin in the center of the granule, and symbol ⊗ represents the dyadic vector product. The related displacements at the infinitely remote boundary follow as with e r , e ϑ , and e ϕ as the base vectors of a spherical coordinate system, also originating from the granule center. Furthermore, the displacement field across the whole domain reads as = (r, ϑ, ϕ) = ξ r (r) sin 2 ϑ cos 2ϕ e r + ξ ϑ (r) sin ϑ cos ϑ cos 2ϕ e ϑ + ξ ϕ (r) sin ϑ sin 2ϕ e ϕ .
The displacement field is then inserted into Equation (11). The resulting strain field enters again the constitutive relations, Equations (12) - (14), and the obtained stress field gives access, via the equilibrium condition, Equations (6) -(9), to a set of ordinary differential equations, which can be solved for the components of the displacement field, see the Supplementary Material. The general solutions of the set of differential equations read with β(α i,j ) defined as The underlying parameters, i.e., α i,j , P i,11 (α i,j ), and P i,12 (α i,j ), are defined in the Supplementary Material. Note that, in line with Bertrand and Hellmich (2009), Equation (33) -(35) are subsequently employed for quantification of the displacement field in the anisotropic constituents only (i.e., solely in the bone phase). For the isotropic constituents within the RVE defined on hierarchical level III, i iso = gran, congl, the displacement field solutions originally given by Hervé and Zaoui (1993) are employed: For determination of parameters Ŵ µ i,j , see the Supplementary material.
For estimating strain states in the macropores, the standard Eshelby matrix-inclusion problem (Eshelby, 1957) is considered, where the infinite matrix is subjected to the same strains as those in the Hervé-Zaoui problem. Under these conditions, the strains in a spherical inclusion with the stiffness properties of the macropores, C macroφ , which is embedded in a matrix with the stiffness C congl of the overall scaffold-bone conglomerate material, read as with S congl sph as the Eshelby tensor of a spherical inclusion embedded in a matrix with stiffness C congl sph , see the Supplementary Material for more details. Notably, the macropores are, as microand mesopores, considered to be drained, i.e., C macroφ = 0. The phase stresses corresponding to the phase strains given through Equation (40)  Homogenization of bulk and shear moduli of the overall conglomerate will be performed on the basis of (i) the microdisplacement, microstress, and microstrain fields in the granule, bone, and macropore phases as given through Equations (23) -(29), (33) -(39), as well as (40) and (41), of (ii) the stress and strain average rules applied to the RVE III of Figure 2, and of (iii) the definition of the bulk and the shear modulus in the terms of macroscopic stress and strain measures. The aforementioned stress and strain average rules, enforcing equilibrium and kinematic compatibility within the RVE III, read as and with the location vector x labeling points inside the RVE III. By definition, the bulk modulus describes the volume change of a material volume if it is subjected to hydrostatic pressure. Accordingly, the bulk modulus of the scaffold-bone compound material with a stiffness according to Equation (15) is defined as where congl,m is the mean macroscopic stress of the scaffold material, congl,m = tr congl /3, and E congl,vol is the macroscopic volumetric strain of the scaffold material, E congl,vol = tr E congl ; "tr" denoting the trace operator. Dividing the RVE III into the material phases introduced in Section 2.1 and Figure 1, the mean conglomerate stresses and the volumetric conglomerate strains can be expressed by means of the stress and strain average rules defined by Equations (42) and (43) and the corresponding mean stress follows as with k macroφ and µ macroφ as the bulk and shear moduli of the macropore phase. The average volumetric phase strains in the granule and bone tissue phases, respectively, are obtained from averaging over strains occurring in the granules, and of the bone tissue coating of the Hervé-Zaoui problem of For the evaluation of Equation (49), it is mandatory that all strain tensors are expressed in the same (Cartesian) coordinate system; also in Equation (49), r i,in and r i,out denote the inner and outer radii of the respective layer: e.g., i = bone gives r bone,in = r gran and r bone,out = r bone , see Figure 2. Analogously, the phase volume averages of the mean stresses in the granule and bone phases follow from the stress field components given by Equations (27)  The mathematical expressions resulting from evaluating Equations (49) and (50) for the granule and the bone phases are given in the Supplementary Material. Finally, inserting Equations (47) -(50) into Equations (45) and (46), and of the resulting expression into Equation (44), yields an implicit equation for calculation of k congl . The macrosopic shear modulus of the scaffold material is governed by another classical relation of continuum mechanics, namely hence µ congl follows from arbitrary components (with i = j) of the deviatoric part of the macroscopic stress tensor, congl,d , defined through congl,d = congl − I congl,m , and the deviatoric part of the macroscopic strain tensor, E congl,d , defined through E congl,d = E congl − IE congl,m . Taking into account that the scaffold material is a composite as defined in Figure 2, with σ d (x) congl gran , σ d (x) congl bone , and σ d (x) congl macroφ as the volume averages of the deviatoric stress tensors across the granular material, bone, and macropore phases, and with ε d (x) congl gran , ε d (x) congl bone , and ε d (x) congl macroφ as the volume averages of the deviatoric strain tensors across the granular material, bone, and macropore phases. The strain tensor components, obtained through insertion of Equations (33) -(35) and Equations (37) -(39), respectively, into Equation (10), are then rotated to a Cartesian base frame, and averaged over the volumes of the individual phases i, Analogously, the volume average of stress tensors across one individual phase, whose components follow from insertion of the phase strain tensors into the constitutive relations, Equations (12) - (14), Given that the here applied loading is purely deviatoric, the deviatoric strain and stress tensors are equal to the respective full tensors, that is ε d ( (56) Due to the purely deviatoric loading, ε d (x) congl macroφ equals ε(x) congl macroφ , reading, when further evaluating Equation (56), as The corresponding average deviatoric stress tensor follows as The macroscopic shear modulus of the scaffold material, µ congl , follows then from insertion of Equations (54), (55), (57), and (58) into Equations (52) and (53), and of the resulting expression into Equation (51). Analogously to the scaffold bulk modulus, µ congl is a function of the scaffold stiffness. This implies that determination of µ congl has to occur in conjunction with determination of k congl , and calls for an implicit scheme.

Animal Studies and Biomaterial Experiments, Revealing Kinetics of Bone Ingrowth and Scaffold Resorption
When applying the hydroxyapatite granule system in vivo, two out of the various phase volume fractions seen in Figure 2 and Section 2.1 significantly evolve over time: (i) the bone tissue volume fraction f congl bone , and (ii) the microporosity φ polyHA micro . Corresponding evolutions can be estimated as follows: As concerns bone ingrowth, we make use of insights gained from X-ray microtomography. Over a time span of roughly 10 weeks after implantation, the increase of the thickness of newly grown bone tissue on implanted biomaterials turned out to be approximately linear, at a rate of approximately k growth = 4 µm/week (Cancedda et al., 2007). Given the near-spherical shape of the granules making up the biomaterial scaffold, the bone volume fraction evolution over time follows as where t is the time variable, t = 0 being the time instant of implantation, and r gran the (average) radius of the granules. In order to define the bone deposition rate that is to be expected for the biomaterial studied in this paper, specific histological animal studies were carried out at the Central Scientific Research Institute of Dentistry and Maxillofacial Surgery (CRID), in Moscow, Russia. These studies were performed in conformity with the relevant institutional guidelines, which are in compliance with respective national laws and policies related to animal care, and which were approved by the Animal Ethics Committee of CRID. Thirty adult Wistar rats (12 weeks old, body weight 250 to 300 grams, and of both sexes) were divided into 2 groups. The rats were allowed free access to food and water ad libitum at all times and were maintained on a 12 h light/dark cycle (lights on from 8 a.m. and 8 p.m.). They were constantly exposed to a temperature of 23 ± 1 • C, and a relative humidity of 60 ± 10 %. All rats were maintained and used in accordance with the guidelines of the Animal Ethics Committee of CRID For surgery, each rat was anesthetized by intraperitoneal injection of Zoletil 50 (Virbac S.A., France) with a dose of 0.2 mL/100 grams of body weight. For the femoral epiphysis model, an incision was made in the skin on the medial side of the thigh and the femoral quadriceps muscle was exposed. The muscle was sectioned longitudinally in its distal third and separated anterolaterally. After exposure of the distal end of the epiphysis of the right femur and peeling of the periosteum, a bone defect was created with a 3 mm hand-held surgical drill under irrigation of NaCl 0.9%. In the test group, carbonate hydroxyapatite (CHA) ceramic granules were inserted into the defect. In the control group, the bone defects were empty and healed under blood clots. Five animals of both groups were sacrificed at each 2, 4, and 8 weeks after implantation. For histological examination, the grafts with surrounding tissues were dissected and fixated in 10% buffered formalin for 24 h. Each graft was sectioned in ten pieces through its midline and embedded in paraffin. Serial 5 µm-sections were deparaffinized, hydrated and stained with hematoxylin and eosin. Photomicrographs of internal sections of each sample were taken by means of a light microscope (Leica DM LB, Germany) and photographed (Sony, Japan). Computer-assisted measurements of the histological parameters were obtained using an automated image analysis system Image-Pro Plus (Media Cybernetics, USA). For each graft at least 5 sections were examined, in each section at least 25 fields of vision were analyzed. The region of interest was taken at 100× magnification and resolution of 2500×1200 pixels, see Figure 3. By measuring the percentage of newly formed bone (area) on the CHA test groups at 2, 4, and 8 weeks after implantation, revealing a bone bone ingrowth rate of k growth = 7±3 µm/week, complying well with the aforementioned previous studies (Cancedda et al., 2007). As concerns hydroxyapatite granule resorption, a recent µCT-based study on tricalcium phosphate biomaterials FIGURE 3 | Histological study of bone regeneration in bone defects optionally containing CHA granules. (A) control group at 2 weeks after implantation-the bone defect is filled with connective tissue; (B) CHA group at 2 weeks after implantation-the granules are surrounded with connective tissue, in some granules osteoid formation can be discerned ; (C) CHA group at 2 weeks after implantation-near bone edges granules are surrounded with woven bone; (D) control group at 4 weeks after implantation-the bone defect is filled with connective tissue and woven bone; (E) CHA group at 4 weeks after implantation-the granules in the center of the defect are mostly surrounded with trabecular and woven bone; (F) CHA group at 4 weeks after implantation (higher magnification)-woven bone is forming around and inside CHA granule; (G) control group at 8 weeks after implantation-newly formed bone can be discerned at the edges of the bone defect; (H) CHA group at 8 weeks after implantation-trabecular bone has formed between CHA granules; (I) CHA group at 8 weeks after implantation-near bone edges granules are integrated in well-formed trabecular bone. (Czenek et al., 2014) revealed that pseudo-physiological conditions favor micropore growth, while pores of larger sizes remain fairly unaffected. The corresponding temporal evolution of the microporosity during granule resorption thus reads as with φ polyHA micro,0 as the microporosity before resorption sets in, and k res as scaffold resorption rate. In order to find numerical values for the resorption rate, the solubility of CHA ceramic granules was studied in a TRIS-HCl buffer solution, exhibiting a pH-value of 7.4 (according to ISO 10993-14-2001) for 21 days at a constant liquid phase volume (thus representing a closed system). The desired pH-value was reached through adding 13.25 g of TRIS (Cat. No: 77-86-1, Sigma-Aldrich) and 125 ml of HCl (Cat. No: 7647-01-0, Aldrich-Aldrich). The solid-toliquid ratio was 0.5 g/100 mL. The development of the calcium concentration over time in the liquid phase was measured using the atomic emission spectrometer Ultima 2 (Jobin-Yvon, France), see Figure 4. The corresponding dissolution rate of the scaffold material can be back-analyzed based on the chemical composition of the material, i.e., Ca 10 (PO 4 ) 6 (OH) 2−2x (CO 3 ) x , with x = 0.05. Additionally, the gained experimental data was fitted by function C Ca 2+ = (0.1073t 2 + 5.18t + 0.3093)/(t 2 + 44.11t + 42.61) (see the dashed graph in Figure 4), in order to get a continuous description of the dissolution progress. Considering the molecular masses of the hydroxyapatite crystals (which constitute the solid part of the granules), we back-calculated, through simple compositional rules, that the dissolution rate amounted to k res ≈ 0.016 w −1 at the start of the dissolution process while, owing to the closed experimental system implying a changing buffer solution in the course of the dissolution process, the dissolution rate decreases thereafter, eventually converging to virtually zero.
FIGURE 4 | Dissolution behavior of CHA granules. Calcium concentration measured in the buffer solution by means of atomic emission spectrometry, reflecting the dissolution behavior of the immersed CHA granule; the cross-shaped markers represent the measured data, whereas the dashed graph represents the fitting of the experimental data by means of function C Ca 2+ = (0.1073t 2 + 5.18t + 0.3093)/(t 2 + 44.11t + 42.61).

RESULTS AND DISCUSSION
3.1. Microstructure-Property Relations: How Micro/Meso/Macroporosity, Crack Density, and Bone Tissue Volume Fraction Govern the Overall Scaffold-Bone Conglomerate Stiffness For elucidating the influence of the scaffold material composition on its stiffness, the latter was micromechanically estimated for the following parameter variations: Furthermore, since stiffness of an isotropic material is often associated with the so-called Young's modulus, we present the computed bulk and shear modulus of the scaffold material in terms of the corresponding Young's modulus, through Accordingly performed stiffness homogenization reveals that the Young's modulus of the scaffold material decreases non-linearly with increasing microporosity, whereas this dependence is the more pronounced the lower the macroporosity, see Figure 5A. Likewise, increasing the mesoporosity leads to a decreasing Young's modulus of the scaffold material, and the decrease is again more significant for lower macroporosities, see Figure 5B.
Comparing Figures 5A,B suggests that the influence of the microporosity on the Young's modulus of the scaffold material is more significant than the influence of the mesoporosityhowever, it should be noted that the technologically relevant range considered in these parameter studies is much larger for the microporosity than for the mesoporosity. An increase of the density of cracks obviously leads to a decreasing stiffness, see Figure 5C. Particularly for low crack densities an increase of the crack density parameter leads to a very steep stiffness loss. For higher crack density parameters, the additional stiffness loss due to further crack densification appears to be far less substantial.
Notably, in the parameter study presented in Figure 5C, crack density parameters up to ǫ = 100 were considered. However, due to the minor dependence of the scaffold material stiffness on ǫ at high crack densities, only the range ǫ = {0, 25} is shown in Figure 5C. Finally, Figure 5D shows the dependence of E congl on the bone tissue volume fraction. Obviously, bone ingrowth leads to a stiffness increase of the scaffold material. Depending on the initial value of the macroporosity, the stiffness increase stops once all of the macropore space is filled with bone tissue. Thus, while the stiffness for φ congl macro = 0.5 is quasi-zero and thus much lower than for φ congl macro = 0.3, the opposite is true once all of the pore space is filled with bone tissue.
Importantly, all dependencies discussed so far are significantly non-linear, and also highly interrelated; e.g., the influence of an increasing crack density differs, in quantitative terms, between an uncracked and a strongly cracked material, and so forth. This implies that for being able to adequately estimate the stiffness of a hierarchical material such as the one studied in this paper, using simplistic empirical relations is certainly not expedient, and may lead to serious misestimations.

Comments on Experimental Validation
The micromechanical representations depicted in Figure 1, and the corresponding homogenization steps, have undergone extensive experimental validation, in the context of various material systems. The self-consistent scheme with needle-shaped solid phases oriented in all space directions and spherical pores, see Section 2.2, has been experimentally validated for various porous hydroxapatite biomaterial systems (Peelen et al., 1978;Akao et al., 1981;de With et al., 1981;Shareef et al., 1993;Arita et al., 1995;Martin and Brown, 1995;Liu, 1998;Charrière et al., 2001;Fritsch et al., 2009); and it has been also corroborated for a wide range of other porous polycrystalline systems, such as gypsum (Sanahuja et al., 2010), or a variety of piezoelectric ceramics, alumina-, circonia-, and silicon-based materials (Fritsch et al., 2013).
The relevance of the Mori-Tanaka estimate for the mesoporous, cracked matrix-inclusion material system (of Section 2.3) can be readily seen from comparison of corresponding results to those of the Finite Element (FE) simulations of Dejaco et al. (2012): µCT scans of the investigated granule reveal a microporosity of φ polyHA micro = 0.445, a mesoporosity of φ gran meso = 0.189, a crack number of N ≈ 3.80 × 10 −6 µm −3 , and an average crack radius of r cr ≈ 260 µm; the latter two yielding a crack density parameter of ǫ avg = 66.85. Inserting this compositional data into the stiffness estimate of Equation (3) yields a homogenized shear modulus of µ hom gran = 78.50 MPa. It needs to be compared to the shear modulus corresponding to the FE-modeled splitting test of Dejaco et al. (2012), which can be retrieved by the analytical formula for a sphere loaded at its poles (Lurje, 1963). The corresponding shear modulus amounts to µ FE gran = 80.63 MPa. The good agreement between µ hom gran and µ FE gran impressively underlines the reliability of the homogenization scheme defined by Equation (3).
The micromechanics of granular assemblies as seen in the top image of Figure 1 are standardly treated by the self-consistent scheme, as has been validated for material systems such as shale (Ortega et al., 2007). This renders the choice of a self-consistent assembly of pores and coated spheres as described in Sections 2.4 and 2.5 as very natural. Of course, an additional direct mechanical testing of bone-scaffold compounds is of general interesthowever, given various technological and ethical challenges, this is clearly beyond the scope of this manuscript. Our philosophy is rather to collect and integrate, by use of latest engineering science developments, the large existing data base in biomaterial science and beyond, in order to bring forth novel design solutions for tissue engineering.

Toward Mathematical Modeling-Based Biomaterial Design
Finally, we present an outlook on how the presented homogenization scheme could be utilized in biomaterial FIGURE 6 | How design parameter combinations affect bone regeneration. Evaluation of bone regeneration kinetics laws, Equations (59) and (60), for r gran = 500 µm, ǫ = 10, k growth = 7 µm/week, and k res = 0.008 week −1 , in terms of microporosity and bone tissues volume fraction evolutions (A), and in terms of the respective Young's modulus (B); model-predicted development of the Young's modulus of the macroscopic scaffold material, E congl , in GPa, during bone regeneration, for varying crack density parameter ǫ = {0, 10, 100}, granule radius r gran = {300 µm, 500 µm, 1000 µm}, macroporosity φ congl macro = {0.3, 0.4, 0.5}, resorption rate k res = {0, 0.008 week −1 , 0.016 week −1 }, and bone ingrowth rate k growth = {4 µm/week, 7 µm/week, 10 µm/week}, depicted for (C) ǫ = 0, (D) ǫ = 10, and (E) ǫ = 100, with the green area ranging from 5 to 19 GPa indicating the targeted stiffness range; graphical representation of successful parameter combinations at (F) t = 5 weeks, and (G) t = 20 weeks, with circle sizes being proportional to the success of specific parameter combinations and line thicknesses being proportional to the success of specific parameter combinations. design. In particular, the important question of how the stiffness of the scaffold material, once implanted, evolves during bone regeneration is addressed-notably, bone regeneration involves growth of new bone tissue, with rate k growth , and resorption of hydroxyapatite crystals, with rate k res , see Section 2.6. First, the development of the material's macroscopic Young's modulus is studied when considering a specific set of parameters, namely, r gran = 500 µm, ǫ = 10, k growth = 7 µm/week, and k res = 0.008 week −1 . Corresponding evaluation of the bone regeneration kinetics laws considered in this work, see Equations (59) and (60), shows that the development of the bone tissue volume fraction depends on the macroporosity, while the resorption-related increase of the microporosity depends solely on the resorption rate, see Figure 6A. The Young's modulus increases due to the addition of new bone tissue in the macroscopic pore space up to the point where all pore space is occupied by bone tissue, after which the Young's modulus slightly decreases due to still progressing resorption of hydroxyapatite, see Figure 6B. Thus, although at t = 0 the macroscopic stiffness increases with decreasing macroporosity, the eventual stiffness after filling up the macroscopic pore space by new bone tissue is actually proportional to the macroporosity.
Extending the elucidation of how the macroscopic stiffness is influenced by the various design parameters of the scaffold material leads actually to a multi-factorial task. Here, we consider the following parameter variations: ǫ = {0, 10, 100}, r gran = {300 µm, 500 µm, 1000 µm}, φ congl macro = {0.3, 0.4, 0.5}, k res = {0, 0.008 week −1 , 0.016 week −1 }, µm/week, 7 µm/week, 10 µm/week}. Then, the arising 243 parameter combinations yield as many corresponding developments of the Young's modulus E congl over time, see Figures 6C-E. Across the human mandible, for which the here studied scaffold material has been developed, a significant variation of the apparent density of the bone organ, and consequently of the stiffness has been observed (Kingsmill and Boyde, 1998;Swasty et al., 2009;Daegling et al., 2011). E.g., in (Daegling et al., 2011), the Young's modulus of mandibular bone is revealed to vary between 5 and 19 GPa-this range of targeted stiffness is indicated in Figures 6C, D by green color, whereas too low or two high stiffness are indicated by red color. It is instructive to evaluate which specific parameter values and which specific parameter combinations lead to the required scaffold stiffness within a reasonable duration. Figures 6F,G show, on the one hand, how often (out of the respectively possible 81 parameter combinations) a specific parameter value has led to a stiffness within the targeted range-the larger (and the darker) the circles related to a specific parameter value, the more often this parameter has led to an appropriate Young's modulus. On the other hand, the lines between two specific parameters indicate how often a specific parameter combination was successful-the thicker (and the darker) the lines, the more often this parameter combination has led to an appropriate Young's modulus. It is striking that after 5 weeks of bone regeneration, the crack density is the by far dominating parameter, i.e., only uncracked materials allow to reach the targeted stiffness, while all other parameters appear to be of minor importance, see Figure 6F. After 20 weeks, however, a low crack density still plays an important role, but also a low granule radius (r gran = 300 µm), a low resorption rate (k res = 0), and a high bone growth rate (k growth = 10 µm/week) turn out as essential for the stiffness development of the scaffold material. In terms of parameter combinations, particularly the combination of ǫ = 0 and k res = 0 turns out as crucial factor for the stiffness development.
Considering that the scaffold resorption rate and the bone ingrowth rate are kind of tunable properties, e.g., by slightly changing the chemical composition of the hydroxyapatite needles, or by adding bone morphogenetic proteins, it is interesting to study the stiffness development for continuous variations of these two quantities. For this purpose, a scaffold material is studied with ǫ = 10, φ congl macro = 0.4 and r gran = 500 µm, whereas k res and k growth are varied in the ranges of k res = [0, 0.016 week −1 ] and k growth = [4 µm/week, 10 µm/week]. Then, our model suggests that after 5 weeks none of the considered combinations of k res and k growth leads to E congl ≥ 5 GPa, that after 10 weeks, a combination of low resorption rate and high ingrowth rate leads to sufficient stiffness, and that after 20 weeks virtually all bone growth rates (within the considered range) lead to sufficient stiffness as long as the bone resorption rate remains at a low level, see Figures 7A-D. Figure 7E shows at which time instants E congl ≥ 5 GPa is reached for specific combination of k res and k growth . Similar computations can be performed for other combinations of ǫ, r gran , and φ congl macro , in order to work out the appropriate parameter space.
Actually utilizing the stiffness prediction tool proposed in this paper for scaffold optimization, requires however more detailed knowledge on the bone ingrowth and scaffold resorption kinetics. In the here presented studies, linear kinetics laws were assumed. In the physiological environment of the implanted scaffold material, long-term occurrence of such linear behavior is however unlikely, for several reasons: Bone ingrowth is a process which typically occurs in degressive fashion, see e.g., (Hing et al., 2004;Cancedda et al., 2007); on the one hand because of purely geometric reasons, namely due to decreasing area on which new bone tissue can be laid down, and on the other hand because of the decreasing availability (with decreasing porosity) of the biological factors necessary for initiating bone apposition. Also, specifically designing a scaffold, for a particular patient, requires information on the location of insertion, e.g., density distributions gained from CT imaging may be sufficient to estimate the site-specific stiffness distribution (Hellmich et al., 2008). Considering these additional influences, namely both the biological environment and the mechanical environment which are to be expected in vivo, is however beyond the scope of this paper.

CONCLUSIONS
In this paper, we have presented a new mathematical model for determining the stiffness of a hydroxyapatite-based granular biomaterial, through upscaling of compositional and morphological information known on three distinct observation scales. The model predictions are, from a qualitative point of view, plausible. The distinctively non-linear dependencies of the derived elastic constants on underlying nano-, micro-, and macro-porosities, as well as on the extent of cracks emerging on the micro-scale underline the necessity of using sophisticated mathematical models when reliable stiffness estimates are required.
The proposed model provides substantiated outlooks as to how specific design parameters influence the macroscopic stiffness of the scaffold-bone compound. In order to reach a physiologically relevant Young's modulus of mandibular bone, typically ranging from 5 to 19 GPa (Daegling et al., 2011), granules with the lowest possible crack density need to be produced, especially at early stages of bone regeneration. Later, other design parameters, such as a low granule radius or a high bone formation rate become relevant as well. In such a way, as demonstrated in Section 3.3, optimal design parameters or even combinations of design parameters can be worked out, paving the way to mathematical modeling-driven optimization of the scaffold-bone compound's performance.
Furthermore, we plan to couple the here presented mechanical model with computed tomography (CT) image-to-mechanical properties conversion techniques, as demonstrated for different materials in (Scheiner et al., 2009;Dejaco et al., 2012). Such conversion techniques allow for direct interpretation of the gray value distributions that constitute CT images in terms of the corresponding distributions of mechanical properties. Finally, once such three-dimensional distributions of the scaffold stiffness and the surrounding bone matrix are computed, the model could be further coupled to recently developed mechanobiology models (Scheiner et al., 2013), eventually allowing to assess the bone regeneration progress for prescribed mechanical loading regimes. These model upgrades will allow to design patientspecific scaffold materials and structures, which is believed to entail a significant increase of the efficiency of such materials.

AUTHOR CONTRIBUTIONS
SS: Conceptual model design, model implementation, performance of all numerical studies, interpretation of the results, as well as drafting and all tasks related to finalizing the manuscript. CH: Conceptual model design, interpretation of the results, and all tasks related to finalizing the manuscript. VK: Processing of the studied material, contributions to interpretation of the experimental results, final proofreading of the manuscript. AG: Implementation of histological and dissolution studies, final proofreading of the manuscript.