Hierarchical Biomechanics: Concepts, Bone as Prominent Example, and Perspectives Beyond

In an interesting, recently published review article, Alber et al. [1] propose the following definition of multiscale modeling: “Multiscale modeling is a successful strategy to integrate multiscale, multiphysics data and uncover mechanisms that explain the emergence of function.” In this context, Alber et al. [1] refer to “small-scale” constitutive equations, and “large-scale” conservation laws as the general physical frame in which the multiscale models, being themselves supported by machine learning if necessary, may be developed. Thereby, the small-scale constitutive equations may be adopted from the classical engineering field, by considering elastic, viscoelastic, or elastoplastic behavior. As concerns living materials, the constitutive equations may be more complex, and concerning the “large-scale laws”, biological applications may require rethinking “of the underlying kinetics, the balance of mass, and the laws of thermodynamics”, and also the inclusion of “the biological, chemical, or electrical fields that act as stimuli of this living response”. We here embrace the aforementioned focus on integration of multiscale/multiphysics data, and we share with Alber et al. [1] a deep root in engineering mechanics, as evidenced through their association of small-scale differential equations with classical mechanical behaviors, such as elasticity [2], viscoelasticity [3,4], and elastoplasticity [5,6]. Accordingly, we set, in this review paper, the primary focus on multiscale mechanics modeling. This choice appears natural, given the importance of mechanics-related aspects in the biomedical field, in particular when focusing on multiscale modeling. However, we wish to complement the account of the state of the art by Alber et al. [1] in multiscale modeling of biomedical materials and systems, by adopting an extended theoretical frame. More specifically, we set our focus not so much on one small (“constitutive”) versus one large (“conservation laws-related”) scale, but on the numerous intermediate scales, where one type of constitutive behavior is effectively transformed into another type of constitutive behavior. At the same time, we wish to review and discuss how the constitutive behavior at the smallest relevant scale of a biomechanical system may be either governed by largely invariant (“universal”) continuum mechanical properties, or even from atomistic mechanics, in both cases avoiding any empirical parameter fitting. It is this focus on the intermediate scales and on the corresponding successive scale transitions, which defines the (bio-)mechanical science-related literature reviewed and discussed in this paper, referred to by the term “hierarchical biomechanics”. Still, we do not restrict ourselves to multiscale methods for pure mechanics, but we also address the exploration of the cross-roads of mechanics/physics and biology, as it is well-known that biologically driven processes may influence mechanical properties, and that, in turn, mechanical loading may influence the biological environment. In this context, the present review is completed by multiscale modeling concepts which, in a theoretically consistent manner, obey both the traditions of mechanics/physics (by employing field theories, interaction/constitutive laws, and geometrical compatibility), and of biology (by employing cell population dynamics, and cell biology models).

As a rule, biomedical materials and biological tissues are hierarchically organized, encompassing multiple pore spaces and other organizational patterns at different length scales [7,8]. Therefore, the ad hoc adoption of the material point concept inherent to traditional solid mechanics is not straightforward, and in most cases even questionable. This suggests prudence and good reasoning whenever using standard continuum mechanics quantities associated with the so-called material point, such as the Cauchy stress tensor (and associated notions of mechanical strength), various types of strain tensors (in particular the Green–Lagrange or the linearized strain tensor), or corresponding stress–strain laws (also referred to as constitutive laws), in the context of highly complex biological and biomedical materials. A good starting point for tackling the corresponding challenges is to resort to well-accepted, theoretically formulated and practically useful concepts which originate from geotechnical engineering, and which actually form the very basis of the scientific field called poromechanics [9]—, i.e., mechanics of porous media, rather than mechanics of pure solids.

A conceptual crystallization point in this context is the so-called representative volume element (RVE), which may be traced back to Biot's famous 1941 treatise on consolidation [10]—the latter is, in Biot's words, a “phenomenon […] whose mechanism is known to be in many cases identical with the process of squeezing water out of an elastic porous medium.” For the appropriate representation of such an elastic porous medium, Biot proposes: “Consider a small cubic element of the consolidating soil, its sides being parallel with the coordinate axes. This element is taken to be large enough compared to the size of the pores so that it may be treated as homogeneous, and at the same time small enough, compared to the scale of the macroscopic phenomena in which we are interested, so that it may be considered as infinitesimal in the mathematical treatment.”

We see that Biot defines a small element, at the same time, both as a point at the macroscopic scale and as a volumetrically extended (here cubic) entity at the microscopic scale. Explicitly assigning two scales to one and the same element goes impressively beyond the material point concept of classical solid mechanics, and hence, is not only qualified as the birth of poromechanics, but, from a more general viewpoint, marks the actual starting point of multiscale modeling in the engineering sciences—long before this term eventually gained ubiquitous popularity, both in the traditional engineering sciences and beyond, especially so in biomedical engineering. The concept of such a simple RVE consisting of a solid matrix and a fluid-filled pore space entails a fundamentally new constitutive property. The latter is called Biot coefficient in the isotropic case [10], and Biot tensor (of second order) in the general anisotropic case [11], and it links pore pressure variations, under constant macroscopic strains, to macroscopic RVE-related stresses. The RVE-concept implies repeated introduction of one and the same type of material property, namely, at different length scales: Elastic properties are assigned to the solid matrix at the microscopic scale and to the overall porous medium, i.e., to the entire RVE. As an additional complication, the overall, RVE-related elastic properties depend again on the action of the fluid, with two important limit cases, as first described by Gassmann [12]: (i) the case of empty pores or of pores with pressures governed from outside of the RVE (drained stiffness), and (ii) the case of a closed RVE with the fluid being trapped inside the RVE (undrained conditions). Corresponding, meanwhile classical relations were applied to soils, rocks, and shales, and nurtured important technological progress in geotechnical and petroleum engineering throughout the next four decades. It was not before the late 1990s, that such basic poromechanical relations entered the bio-engineering realm. This was largely the merit of Cowin and coworkers, who dealt with the poroelasticity of bone [13,14], and transferred, very successfully, the consolidation problem into biomechanics, by discussing the role of fluid flow in stimulating biological cells residing in pores, one of the key topics of this review.

Independently of the early poromechanical developments, the idea of an RVE was also embraced by the community dealing with reinforced solids. Going beyond the poromechanicians' conceptualization of a fluid pressure at the microscale, the reinforced solids community introduced the entire stress and strain fields at the microscale, thereby establishing, in the 1960s, a new scientific field called composite mechanics or continuum micromechanics [15]. The key concept of the latter is that a macroscopic material point is associated with a full continuum mechanics boundary value problem, formulated on an RVE defined at the microscopic scale. In terms of length scales, this RVE coincides with that of poromechanics. More specifically, the size of the pores is now taken by the size of any type of (solid) inhomogeneities, and the boundary of the RVE is subjected to a microscopic displacement field which is governed by homogeneous macroscopic strains. These boundary conditions readily imply a strain average rule [16,17]: The spatial average of kinematically compatible microscopic strains is equal to the macroscopic strain. Alternatively, if the boundary of the RVE is subjected to microscopic traction forces which are associated, via Cauchy's fundamental theorem, with one and the same macroscopic stress tensor, the following stress average rule applies: The spatial average of equilibrated microscopic stresses is equal to the macroscopic stress. Accordingly, one of the pioneers in the field, Rodney Hill, writes in 1963, when explaining what he calls “the representative volume” [16]: “This phrase will be used when referring to a sample that (a) is structurally entirely typical of the whole mixture on average, and (b) contains a sufficient number of inclusions for the apparent overall moduli to be effectively independent of the surface values of traction and displacement, so long as these values are macroscopically uniform. That is, they fluctuate about a mean with a wavelength small compared with the dimensions of the sample, and the effects of such fluctuations become insignificant within a few wavelengths of the surface. The contribution of this surface layer to any average can be made negligible by taking the sample large enough.” Hill has also complemented the stress and strain average rules by a third average rule associated with (rates of) work, standardly referred to as Hill's lemma: The average (rate of) work done by equilibrated microstresses on kinematically compatible (rates of) microstrains is identical to the (rate of) work done by the macroscopic stresses on the macroscopic strains.

In this context, the following practical question arises: Which ratio between the characteristic length of the RVE and that of the typical inhomogeneity within the RVE is needed to keep the aforementioned fluctuations sufficiently small, i.e., to make homogenization over microheterogeneity admissible and useful? This question had kept the scientific community busy for some time, until Drugan and Willis [18] provided a particularly satisfying answer, by considering the ergodicity principle and employing ensemble averaging. More precisely, Drugan and Willis [18] consider the RVE as an ergodic system where “local configurations [of microscopic quantities, such as polarization stresses quantifying the fluctuation of elasticity tensors] occur over any one specimen with the frequency with which they occur over a single neighborhood in an ensemble of specimens.” Hence, the spatial average of microstrains over one material sample representing the small material element in the sense of Biot and Hill becomes equivalent to the average of the microstrain at one specific microscopic point over an ensemble of very many samples representing the same type of material. Subsequent comparison of the mathematical solution for the stress–strain relations of the homogenized single sample with that of the ensemble averages reveals a surprising result: An RVE exceeding only twice the diameter of spherical reinforcements defining the microheterogeneity size already allows for an accurate prediction of the overall elasticity, as quantified by an error margin of only 5%. Only one year later, this somewhat surprising result was heuristically confirmed by Monte Carlo computations on cubic unit cells hosting very many different disordered arrangements of spheres subjected to periodic boundary conditions [19], and the same scale separation factor holds for parallel, infinitely long cylindrical inhomogeneities [20]. The rather small RVE size also indicates that the boundary of an RVE undergoing homogeneous (i.e., Hill-Hashin) boundary conditions must be carefully chosen: It needs to circumvent local environments with high elasticity fluctuations, as explicitly noted by Dormieux et al. [21] when discussing pores as “zero-stiffness reinforcements”.

The aforementioned compact RVE sizes impressively underline the role of continuum (micro-)mechanics as a particularly efficient, robust, and versatile tool for materials with pronounced microscopic organizations, such as biological tissues. However, similar to the situation with poromechanical theories, RVE approaches hardly entered the biomedical field before the 1990s, and it was not before the 2000s until they started to flourish. They did particularly well for the intricate hierarchical organization of bone tissue where classical anatomical terms could be well assigned to RVEs of different sizes: trabecular bone to an RVE with a size of several millimeters [22–30], cortical bone to an RVE with a size of several hundred micrometers [20,23,31–46], extravascular bone material to an RVE with a size of hundred micrometers [35,46,47], extracellular material to an RVE with a size of tens of micrometers [32,35,39,41,47–51], and finally both mineralized collagen fibrils [34,35,47,48,52–54] and extrafibrillar polycrystals [34,35,47, 52,55] were assigned to RVEs with a size of about hundred nanometers. The latter two RVE types contain the elementary building blocks (or “major components” [7]) of bone: i.e., nanometer-sized hydroxyapatite crystals and crosslinked type I collagen molecules [52]. Notably, when using semi-analytical approaches, as discussed in more detail in Sec. 3, several of the aforementioned RVEs were straightforwardly combined to multistep homogenization schemes, sometimes referred to as “Russian doll models” [23,24,32,34,35,39,41,48–50,56–67]. In this case, the microheterogeneities within an RVE are represented by yet smaller RVEs, with the size of the latter being smaller or equal to the aforementioned microheterogeneity size [47]. An illustrative example is depicted in Fig. 1, featuring the hierarchical organization of cortical bone, from the macroscopic, structural scale down to the scale of single collagen molecules, hydroxyapatite crystals, and fluid-filled pore spaces, together with the corresponding multistep micromechanical representation. Again, the actual realization of homogeneous stress or strain boundary conditions imposed onto a sufficiently RVE may practically imply carefully chosen, somewhat wavy RVE boundaries deviating from Biot's originally envisioned element (being of cubic shape), as can be seen from Level 5 shown in Fig. 1.

While bone appears as the most intensively studied biological tissue, several other biomedical materials have been studied by RVE approaches as well, including arterial wall tissue [73,74], muscle tissue [75], skin tissue [76–78], heart tissue [79], dentin [80–82], lung tissue [83,84], or synthetic tissue engineering materials [85–90].

The general nature and usefulness of the RVE being set now, it is time to continue with connecting microscopic to macroscopic physical quantities i.e., with scale transition: upscaling and downscaling.

The quest for efficient and/or accurate solutions for microscopic boundary value problems associated with RVEs has stirred enormous scientific efforts and led to many applications of standard computational mechanics-based approximation tools, such as the Finite Element (FE) method [20,29,76,83,85–94]. In the latter context, periodic boundary conditions are usually prescribed in a detailed representation of the RVE's microstructure. The latter may be (i) geometrically well-defined, as is the case for synthetic tissue engineering scaffolds [85–87,89,90]; (ii) generated through suitable algorithms imitating key morphological features, such as fractal trees in skin microvasculature [76] or hexagonal staggered prisms in dental enamel [91]; or (iii) reconstructed from imaging data provided by techniques such as scanning acoustic microscopy [20,95] or microcomputed tomography (micro-CT) [29,83,88,92–94,96]. While appealing from a “ready-to-use” perspective which does not require deeper conceptual scrutiny, such approaches require a “full” description of the microstructure. Given the corresponding high computational demands, the latter are typically restricted to single-step homogenization procedures. Hence, the material properties entering the FE models on the integration point level either need to be known or suitably guessed, or they follow from additional, then analytical or semi-analytical homogenization and scale transition techniques “digging” deeper into the micro-and nanostructures of the investigated material systems [66,88,96].

With the FE method, the CPU time scales with the square of the degrees-of-freedom, which may render this method inefficient for very large complex problems, where, in addition, meshing may become a challenge on its own. As a remedy, both the scaling and meshing problems can be elegantly overcome by numerical methods based on the Fast Fourier Transform (FFT). The latter methods rely on the smart use of Green's functions [97–99], i.e., of the solution for the displacement field throughout an infinitely extended homogeneous elastic body subjected to a singular point load of unit intensity. First, the microheterogeneous linear elasticity field is equivalently represented as the superposition of a homogeneous elastic field and a field of microheterogeneous eigenstresses, so-called polarization stresses. Setting the microstress and microstrain fields of both representations equal to each other delivers the polarization stresses as the tensor product of the microheterogeneous deviations of the micro-elasticity field from the properties of the homogeneous elastic body on the one hand, and the microstrains on the other hand [100–102]. This opens the way to incorporating solution schemes for the Lippmann-Schwinger equation—which was originally proposed in the framework of quantum mechanics [103]—into the conceptual framework of continuum micromechanics, particularly in the context of FFT analyses. The latter was pioneered by Moulinec and Suquet [104] and refined thereafter by Brisard and Dormieux [105]. Only very recently, namely, since 2018, these elegant solutions have made their way into the image-based micromechanics of cortical bone [38,67]. Such approaches outperform FE-based homogenization approaches by orders of magnitude in terms of computational efficiency—however, they still rely on “complete” microstructural information. This might be the reason why yet another approach, namely, that based on so-called material phases and estimates for homogenized mechanical properties, has gained much more popularity in the theoretical and applied mechanics field; and substantially propelled true multiscale approaches in biomechanics and bio-engineering.

The premise for phase-based estimation of RVE-related properties is the explicit awareness that the resolution of the microstructure in infinite completeness is eventually impossible. Accordingly, this resolution is not even aimed at, but the focus is a priori set at microstructural features potentially governing the homogenized mechanical behavior of the RVE. To that end, homogeneous subdomains, referred to as material phases, are introduced within the RVE [15,17,106]. Each material phase is characterized by homogeneous mechanical properties (such as a phase elasticity tensor, instead of a micro-elasticity field), by a characteristic shape, and by the volume fraction it occupies within the RVE.

As the first key modeling step in phase-based homogenization, each phase is approximated by a homogeneous ellipsoidal elastic inhomogeneity, and each of those inclusions is embedded into the same type of infinitely extended auxiliary matrix subjected to homogeneous strains at its infinitely remote boundary, see Fig. 2. This auxiliary matrix exhibits elastic properties which are different from those of the inhomogeneity. The interest in such auxiliary matrix-inhomogeneity problems arises from the existence of corresponding analytical and semi-analytical solutions. They were provided by Eshelby [107] and by Laws [108], for inhomogeneities embedded in isotropic and anisotropic matrices, respectively. As was the case with the FFT solutions described in the last paragraph, also the solutions for Eshelby's matrix-inhomogeneity problems were derived by combining the polarization stress concept with the method of Green's functions. These solutions provide the remarkable result that the strains in the inhomogeneity are actually homogeneous and depend on the inhomogeneity shape, the elastic stiffness contrast between the inhomogeneity and the surrounding matrix, as well as on the remote matrix strains. If the investigated RVE is predominantly made up by a matrix phase in which inclusion phases are dilutely embedded, the RVE-related macroscopic strains can be approximated by the auxiliary matrix, and the macro-to-microstrain transition is already completed. Insertion of the phase strains into the phase-specific linear elastic laws yields phase stresses and averaging the latter over the RVE eventually gives access to macroscopic stresses; leading to the homogenized elasticity tensor associated with the so-called dilute homogenization scheme. Modeling of “dilute suspensions” can be traced back to Einstein's 1906 treatment of rigid spherical inclusions in a Newtonian viscous fluid [109], and already appears as a broad review topic in Ref. [106].

However, for the frequently encountered cases where the individual phases are equally dominant in filling the RVE and consequently interact with each other, a second key modeling step is necessary, as sketched by Mori and Tanaka in 1973 [110] and rigorously formulated by Benveniste in 1987 [111]: Namely, the matrix-inhomogeneity problems-related phase strains are inserted into the strain average rule, thereby providing a link between the remote auxiliary matrix strains and the macroscopic strains associated with the RVE, see Fig. 2. This link is then combined with the aforementioned solutions for the matrix-inhomogeneity problems, yielding relations between phase strains and RVE-related macroscopic strains. The latter relations are multilinear and hence expressed by so-called strain concentration or downscaling tensors of fourth order. Use of the phase-to-RVE strain relations in the phase-specific microscopic linear elastic law, and subsequent averaging of corresponding microstresses over the RVE yields the macroscopic (homogenized) RVE-related elasticity tensor, as an analytical or semi-analytical tensor function involving the elastic properties of the phases, their shapes, their volume fractions, and the stiffness of the auxiliary matrix [15], see also Fig. 3. Choice of the latter completes the elasticity upscaling problem, with two main options: Choosing the matrix as an individual phase with given stiffness and volume fraction leads to the so-called Mori-Tanaka scheme [110], which is appropriate for matrix-inclusion composites. The second option is to assign the homogenized elastic properties and a zero volume fraction to the matrix. This leads to an implicit equation for the homogenized elastic properties. Its iterative solution is called self-consistent scheme. This scheme is appropriate for polycrystalline materials where all phases are in intimate mutual contact. Self-consistent schemes in the context of continuum mechanics can be traced back to the 1950s and 1960s [112–114]. However, they effectively entered bio-engineering only in the early 2000s, in the context of modeling the porous polycrystals (or mineral foams) [115] and the mineralized collagen fibrils [49] present in the ultrastructure of bone. From that time on, also the Mori-Tanaka homogenization scheme was successfully applied to various levels across the hierarchical organization of bone, namely, for the representation of wet collagen (Level 1 in Fig. 1), of the extracellular bone matrix (Level 3 in Fig. 1), of the extravascular bone matrix (Level 5 in Fig. 1), and of the macroscopic bone material (Level 6 in Fig. 1). The corresponding semi-analytical tensorial expressions are given in Fig. 3.

Bio-engineering applications were a major driving force for developing Eshelby problem-based homogenization schemes for material properties beyond the elastic regime. The corresponding classical extension is that from elasticity to viscoelasticity, i.e., to assigning Boltzmann-type linear viscoelastic behavior to the material phases which per se remain fixed throughout time. In the Laplace-Carson space, such a linear viscoelastic behavior appears as a succession of formally linear elastic material behaviors, for each of which the classical Eshelby problem-based homogenization schemes described in the two preceding paragraphs can be employed, as explained in detail by Laws and McLaughlin in 1977 [116]. Effective application of this type of theory in bio-engineering took place only in 2014 [117], by assigning viscoelastic properties to the hydrated extrafibrillar mineral crystals in the ultrastructure of bone, in order to explain hydration- and mineralization-dependent creep of macroscopic bone material. In this context, the computational endeavors really start when back-transforming the homogenized properties into the time domain; e.g, by means of the Gaver-Wynn-Rho algorithm [118].

Another straightforward extension beyond elastic homogenization concerned perfectly brittle materials where failure at the microscopic level is directly associated with overall macroscopic failure of the RVE. Accordingly, the insertion of the phase stress-to-RVE stress relation (expressed by fourth-order stress concentration tensors) into a microscopic failure criterion directly yields a macroscopic failure criterion, which considers not only the strength properties of the weakest phase, but also the mechanical interaction of phases, driven by their shape, arrangement, and elastic stiffness contrast. However, when comparing the evaluation of corresponding mathematical formulations with experiments (see Sec. 4), the representation of the RVE's microstructure in terms of a finite number of phases generally turns out as a too crude approximation. As a remedy, the RVE may be represented by infinitely many nonspherical phases being oriented in all space directions, as indicated in Fig. 1, Level 2B, and labeled by polar and azimuthal angles [119]. Corresponding two-dimensional integrals over the unit sphere are then solved by appropriate approximation schemes such as Stroud's formulae [120]. This type of multiscale modeling has provided unprecedented insight into the strength characteristics of man-made hydroxyapatite biomaterials [121–123]—and more generally, of a large number of ceramic material systems [124,125].

However, what really broadens the application range of Eshelby-problem-based homogenization, is the explicit introduction of eigenstresses and eigenstrains in the true physical sense, i.e., the introduction of strains and/or stresses which are fully independent of elastic, i.e., nondissipative, phenomena. In the context of biomedical materials, two types of such eigen-quantities are of particular interest: plastic strains associated with ductile behavior, and pore pressures arising from fluids entering and leaving the various pore spaces found throughout the materials' hierarchical organizations. Important averaging rules for eigenstresses and eigenstrains go back to Levin [126] and follow from repeated use of Hill's lemma, as explained in Refs. [15,127]: The macroscopic eigenstress tensor is the average over the RVE, of the tensor contraction between the microscopic eigenstress tensor and the strain concentration tensor; and the macroscopic eigenstrain tensor is the average over the RVE, of the tensor contraction of the microscopic eigenstrain tensor and the stress concentration tensor. The aforementioned eigenstress averaging rule implies a most interesting link between the micromechanics of media with random microstructures and the classical poroelasticity theory of Biot, as lined out in Ref. [21]: Considering a two-phase composite consisting of an elastic matrix and a pore space with eigenstresses of the format pore pressure times the second-order unity tensor yields the Biot tensor as the pore space-related strain concentration tensor times the porosity. This has a very important practical implication for the investigation of porous materials: The Biot tensors quantifying the interaction between pore pressure and macroscopic stress under a state free of macroscopic deformation need not be determined from complex poromechanical tests with independent control of pore pressure and macroscopic stress or strain, but can be readily determined once an elasticity homogenization scheme, as the ones discussed further above, has been established. This opportunity has been repeatedly taken in the context of the mathematical modeling of bone, in particular so for the determination of pore pressures in the intermolecular, intercrystalline, lacunar, and vascular pore spaces under various physiologically relevant loading scenarios [46,128–130]. The vascular and lacunar pore spaces host the biological cells involved in bone remodeling: osteoclasts, osteoblasts, and osteocytes, see Sec. 5 and Fig. 9 for more details.

When the interest goes beyond “simple” eigenstress or eigenstrain averaging of a two-phase composite with only one phase being eigenstressed, interesting scale transition solutions concerning eigenstrains and eigenstresses can be obtained from the extension of the classical Eshelby-Laws matrix-inhomogeneity problems to eigenstressed matrix-inhomogeneity problems [15]. The solutions to such problems are again based on Green's functions, now in combination with a homogeneous elastic medium and two types of eigenstresses: the classical polarization stress explained further above, and an additional, physically motivated eigenstress, such as a pore pressure or an eigenstress associated with plastic strains. Benveniste's micro-to-macroscale transitions can then be extended to eigenstressed multiphase media, as outlined by Pichler and Hellmich [131]: Each phase is represented by an inhomogeneity in an infinite matrix subjected to remote auxiliary strains; however, now both the inhomogeneity and the matrix exhibit eigenstresses. Insertion of corresponding phase-specific homogeneous strains into the strain average rule yields a relation between the remote auxiliary strains, the RVE-related macroscopic strains, and the eigenstresses in the inhomogeneity and in the matrix, respectively. In order to eventually relate the eigenstresses in the auxiliary matrix to the RVE-related macroscopic eigenstresses, additional use of Levin's theorem is made, finally yielding the so-called concentration-influence relations. The latter is the cornerstone of Dvorak's “transformation field analysis” introduced in the early 1990s [132–134]: The strains in each phase are not only related to the RVE-related overall macroscopic strains (by the already discussed strain concentration tensor), but also to the eigenstrains in all other phases and in the considered phase itself (by means of so-called fourth-order influence tensors). Homogenization of eigenstressed media was particularly helpful for explaining the strength of bone, arising from ductile sliding between hydrated crystals followed by brittle failure of molecular collagen [35]: Therefore, elementary plastic strains fulfilling the flow and consistency rules of classical plasticity [135] were introduced at the level of the hydroxyapatite crystals at Level 2B of Fig. 1, and thereafter upscaled up to the level of macroscopic bone material. More recent studies based on the representation of Level 2B of Fig. 1 have revealed the elastoplastic failure behavior of the cement lines surrounding secondary osteons in Haversian bone from the surrounding bone matrix [55].

Yet further extension of the use of Eshelby problems for scale transitions in complex materials concerns large deformations, the typical biological example being soft tissues. As compared to hard tissues, this challenge has only recently been tackled. First success approaches dealing with arterial and tendinous tissues [74,136] were mainly based on (i) re-formulation of the Eshelby problem for velocity gradients (rather than linearized strains), allowing for the derivation of concentration relations between macroscopic strain rates on the one hand, and microscopic strain and spin rates, on the other hand; as well as on (ii) thermodynamically consistent hypoelastic material laws [137,138], linking objective (i.e., observer-independent) phase properties. Such Eshelby problem-based approaches may well help to overcome persisting difficulties in hierarchical soft tissue modeling. The latter include (i) the incomplete description of load-induced (micro-)fiber kinematics in the context of deformation gradient-based (hyperelastic) models, requiring rather complex, computationally expensive remedies [139,140]; and (ii) the similarly high computational expense of fiber network models [141–144], which, in addition, do not provide explicit access to the matrices found in between the fibers.

As a rule, the usefulness of all the aforementioned applications of scale transition techniques to biological and biomedical materials needs to be shown through careful and extensive experimental validation. This is the topic of the next section.

According to Popper's famous reasoning on the logic of scientific discovery [145,146], a hypothesis which the human mind comes up with needs to survive as many statistically and physically different falsification tests as possible, in order to eventually reach, after extensive discussions and debates between the peers within a scientific community, the status of a reliable scientific theory. In this spirit, various multiscale micromechanics approaches, as described in Sec. 3, have been carefully validated by means of biomechanical, biophysical, and biochemical experiments [23,35,47,49,52,55,147]. Typically, such experiments concern

• “universal”, i.e. tissue-invariant, mechanical properties of the elementary components making up any tissue belonging to a particular class of tissues (e.g., in the case of bone tissues: molecular collagen type I (making up around 90% of all organic matter found in bone [148]), hydroxyapatite, and water with non-collagenous organic, see Secs. 4B4.2–4.4 for details);

• tissue-specific mechanical properties at any hierarchical level above that of the elementary constituents, tested across the vertebrate animal kingdom, see Secs. 4.2 and 4.4 for details, and

• tissue-specific compositional information which can be converted into phase volume fractions, phase shapes, and phase interaction patterns at any hierarchical level above that of the elementary constituents, again tested across the vertebrate animal kingdom, see Sec. 4.5 for further details.

Based on these tests, the validation strategy is as follows:

• the “universal” properties of the elementary components are associated with the phases which are not anymore resolved into RVEs themselves (i.e. in the case of bone illustrated in Fig. 1 to crosslinked molecular collagen, to hydroxyapatite, and to the intermolecular, intercrystalline, lacunar, and vascular pore spaces);

• the micromechanical model is fed with phase-specific volume fractions arising from biochemical and biophysical tests as detailed further below; and

• the respective model predictions in terms of the mechanical properties associated with the different RVEs are compared to corresponding experimental data on the very same mechanical properties.

This validation strategy has been particularly successful in the context of bone and mineralized tissues; and this concerns different types of mechanical properties, associated with elasticity [23,41,48,49,51,53,61,65,66,115,147], viscoelasticity [117], poroelasticity [128–130], and elastoplasticity/strength [35,55,149]. We note in passing that scattering in material properties typically decreases upon elastic homogenization from smaller to larger scales [150]; and we continue with a more detailed account of the various experimental techniques that have provided data for the aforementioned model validation strategies. The term “universal” in the above paragraphs makes reference to the concept of architectural constraints governing the universal patterns which are omnipresent in the living world, classically discussed in the context of evolutionary biology [151,152]; later, this discussion has been extended toward a mathematical setting for hierarchical biomaterials mechanics [47].

Ultrasonic measurements refer to a very wide range of testing devices, including elastography which “aims at quantitatively image the Young's E modulus, the physical parameter corresponding to the stiffness”, typically throughout larger portions of the organ [153].

Our present focus is different, and concerns the characterization of a particular material sample rather than a larger portion of an organ, and more precisely, the determination of several components of the elasticity tensor defining tissue properties at different length scales. In elastic media, the product of the squared pulse velocity with the sample's mass density gives access to the components of the elasticity tensor [154]. However, ultrasonic tests may also give access to the elastic properties of media which are not behaving purely elastically. In more detail, for impulse frequencies associated with oscillation periods which are much shorter than the characteristic creep times of elasto-viscous materials, the wave velocity is linked to the instantaneous elastic properties of the material, even if the latter may show also inelastic deformations [154]. It is important to note that the actually tested RVE is determined through the wavelength, with the characteristic length of the RVE being five to ten times smaller than the wavelength [155], see Fig. 4. Ultrasonic tests have been very successfully employed for biological and biomedical material characterization, both at the fundamental level of the elementary constituents of a biological tissue class and at various levels of hierarchical organization. In the case of bone, ultrasonic tests at different frequencies give access to the majority of the RVEs depicted in Fig. 1 [23,47,147]: 10 MHz tests [70,156–158] refer to the extracellular bone material (Level 3 of Fig. 1); 2.25 MHz tests [159–161] refer to the extravascular bone matrix (Level 4 of Fig. 1), and 50 kHz tests [162–164] refer to the macroscopic (here trabecular) bone material (Level 5 of Fig. 1).

Ultrasonic tests have also provided access to the elastic properties of the aforementioned elementary components of bone: The bulk modulus of water amounts to 2.3 GPa [165]; and compared to this value, the shear modulus is negligibly small; tests on dense hydroxyapatite powder reveal an isotropic elastic Young's modulus of 114 GPa and a Poisson's ratio of 0.27 [166]. The elementary component “type I collagen” was elastically characterized by Cusack and Miller [167], by an interesting variant of ultrasound transmission concerning very high frequencies, called Brillouin light scattering. Namely, the speed of thermally excited elastic waves is measured in terms of the inelastic scattering of light which these waves emit. The employed 10 GHz signals refer to wavelengths of 300–400 nm, hence characterizing an RVE of some tens of nanometers characteristic length, i.e., to Level 1 in Fig. 1. In more detail, Cusack and Miller tested rat tail tendon, which, upon drying, collapses to a material consisting almost exclusively of type I collagen [168,169] and hence, represents “molecular collagen” phase of the RVE of Level 1 in Fig. 1 [47,49]. This reasoning was further refined [168] by considering xylene imbibition tests [170], showing that dried collagen still exhibits 12% intermolecular porosity in the sense of Level 1 of Fig. 1; and most interestingly, hydration-dependent wave velocity changes reported by Cusack and Lees [171] could indeed be traced back to composite material behavior in the sense of the RVE of Level 1 in Fig. 1. However, the issue might get trickier once a critical water-to-organic mass ratio of about 0.45 is exceeded: then, conformational changes of the molecular collagen itself are probable to lower its intrinsic elastic stiffness.

Totally different from high-frequency elastic waves, the geometrical configurations of biomacromolecules, as observed under the electron microscope, may also give access to their elastic properties. The respective conceptual and theoretical foundations are as follows [172,173], see also Fig. 5: The investigated biomacromolecule is considered as a long, wound thread approximated by a two-dimensional curve of arbitrary (differentiable) shape. Starting from a chosen origin, all the points on the curve are labeled by a contour length variable, and each of these points is assigned to the cosine of the angle between the tangent to the curve at the specific point, and the tangent to the curve at the origin. This process is repeated for many different origins and averaged over all the associated cosine-length relations: As a result, the average over the cosines of angles between any two tangents with a particular contour length between, decays exponentially when considered as function of the aforementioned (increasing) contour length; starting with one, and eventually tending to zero. The length characteristic quantifying this exponential decay is called persistence length. For small angles, the cosine can be expressed by means of the square of the angle, with the latter becoming a function of the contour length and the persistence length, see Fig. 5. The key to relating this statistical distribution of “bending angles” between tangents, to mechanical properties lies in the consideration of different bending energy states, each of which can also be expressed in terms of a bending angle and a contour length. Moreover, the bending energy states follow a Boltzmann distribution governed by the absolute temperature. Equating the Boltzmann distribution-related average of the square of the bending angles, to the corresponding average obtained from the geometrical evaluation of the biomacromolecule's shape, finally provides the bending stiffness (also referred to as flexural rigidity) as a function of the persistence length, the Boltzmann constant and the absolute temperature. When considering, in addition, the biomacromolecule to have a circular cross section, standard beam theory can be used for relating the bending stiffness to the Young's (elastic) modulus of the “material” making up the biomacromolecule. In the case of collagen type I molecule with side chains in a biologically relevant conformation [174–177] a cross-sectional diameter of 1.4 nm, together with a persistence length of 57 nm, results in an elastic modulus of approximately 3 GPa.

The epitome of biomechanical testing is uni-axial quasi-static experiments, where a cuboidal or cylindrical sample is subjected uni-axially to tension or to compression forces, see Fig. 6. Care needs to be taken in order to realize, as far as possible, homogeneous uni-axial mechanical stress states, with a magnitude amounting to that of the applied force divided by the area over which the force is transferred onto the material sample. These mechanical stresses are related to length changes of the sample. The latter may be recorded by very elaborate techniques, such as wide angle X-ray scattering in order to detect length changes at the level line elements which are as small as 0.29 nm [178]. Corresponding tension tests on bovine achilles tendon confirmed the 3 GPa elastic modulus of collagen, as referred to above in Sec. 4.3. Collagen is known to behave, in a fairly good approximation, elastically upon loading, with an abrupt transition from elastic deformations to brittle fracture, as can be inferred from tensile tests on collagen fibers [179]. However, this is not the case any more with more complex composite-type materials such as bone, where elasto-plastic behavior is encountered already at low loading levels [35,55,180].

In order to retrieve elastic properties from such tests, the thermodynamic nature of elasticity needs to be carefully considered [2,3,9,181]: Elastic energy is that portion of the internal energy stored in the material, which can be back-transformed into efficient mechanical work. This transformation happens once a sample is unloaded. Hence, slopes in the loading regime of mechanical tests encompass both (visco-)elastic and plastic portions, and only the unloading regimes guarantee the absence of plastic deformations. Accordingly, moduli derived from the loading portions of stress–strain curves qualify as “(visco-)elasto-plastic”, and only the unloading branches of stress–strain curves are associated with (visco-)elasticity only [182–185]. In this context, we also note that the mineral and water contents are key drivers of bone viscoelasticity: This underlines the importance of the hydrated extrafibrillar polycrystal for the dissipative mechanical behavior of bone [117,182,186,187].

When it comes to strength, there is virtually no alternative to mechanical tests. However, in the context of strength, the load is not lowered so as to get access to (visco-)elastic properties, but it is increased up to mechanical failure of the sample. Accordingly, the strength properties of the elementary components of bone (collagen and hydroxyapatite) were inferred from uni-axial mechanical tests where rat tail tendons [179] and samples of densely compacted hydroxyapatite powder [189,190] have been loaded up to failure. The strength of the extracellular bone material can be gained from micropillar testing [188,191,192], the strength of polycrystalline material making up the cement lines around osteons can be accessed through osteon push-out tests [193], and the strength of macroscopic bone materials has been amply tested over the last five decades (at least). In this context, it is advisable to prefer tests where homogeneous stress states prevail in the sample, such as uni-axial tension tests and uni-axial compression tests [194–211]. For example, bending tests deliver much less reliable results as major hypotheses concerning material and structural behavior need to be made. Normally, such hypotheses have been readily borrowed from engineering steel construction—however, as actually evident from abundant mechanical and chemical features, bone is obviously not steel.

In this context, demineralization, decollagenization, and/or dehydration tests, in combination with weighing tests, Archimedes' principle, inductively coupled optical emission spectroscopy (ICP-OES), and microsocpic investigations [69,212–218] are the key experimental access to the chemical composition of bone ultrastructure. Surprisingly, across many different bone tissues of differently aged mammalian and avian species, there emerges a clear bilinear relationship between the apparent mass densities of mineral and organic matter (i.e., of the masses of mineral and organic matter divided by the volume of extracellular matrix within which they reside), see Fig. 7 [157,219–225].

When additionally considering X-ray and neutron diffraction tests giving access to the distance between collagen molecules, several “universal” rules concerning the composition and microstructure of mineralized collagenous tissues (such as bone) can be derived: (i) upon hydration, the extrafibrillar space grows proportionally to the fibrillar swelling in nonmineralized collagenous tissues [168]; (ii) within the bone ultrastructure, the average extrafibrillar, and the average extracollageneous mineral concentration are the same [226]; (iii) the mineral and collagen volume fractions within the extracellular bone matrix are linearly related to each other [227]; and (iv) tissue shrinkage upon mineralization can be explained from precipitation of mineral from a liquid solution under closed thermodynamic conditions [228]. These rules yield the tissue-specific volume fractions within the RVEs of Levels 1 to 3 in Figs. 1 and 3; while the pore volume fractions at the higher levels are accessible from light microscopy [225] or computed tomography (CT) [229]. On top of that, the RVE representation of the bone ultrastructure has independently been confirmed by an electrodynamics model based on Maxwell's equations [230].

Experimentally validated material properties retrieved from hierarchical biomechanics models for level 5 of Fig. 1 can also be used for feeding structural mechanics models of organs, such as an FE model of a lumbar vertebra as seen in Fig. 8. Such models are typically derived from CT [231–233]. Thereby, the CT data gives access to (i) the topology of the organ and (ii), via the average rule for X-ray attenuation coefficients [234,235], to heterogeneous fields of vascular porosities throughout the organ. When additionally accounting for the invariance over space and time of the organ-specific chemical composition of the bone ultrastructure when averaged over millimeter-sized domains [225,236–239], and of the correspondingly invariant volume fractions inside the RVEs of the ultrastructure level and below, the aforementioned homogenization schemes provide inhomogeneous mechanical property distributions throughout the FE models. Eventually, this allows for the estimation of patient- and activity-specific safety factors of a load-bearing organ; elucidating that mild physiological loading of a lumbar vertebra is associated with a safety factor of around five [231], see Fig. 8. This value is fully consistent with the ultimate forces obtained from mechanical tests on whole vertebrae [240–243] as well as with loadings encountered in the context of extreme sports activities [244]. Similar types of structural modeling have been realized for various other systems, including human and murine femora [233,245,246], human mandibles [232,247], or different tissue engineering constructs [88,96,122,123,248].

As is evident from the preceding sections, hierarchical biomechanics has most maturely evolved in the case of bone tissue, with first interesting approaches appearing for other tissues as well. All these results concern mechanical properties, and how they are related to tissue composition and organization. However, hierarchical modeling may still enter the picture in a yet broader scope, namely, when studying the interplay between mechanics and biology, either at the hierarchies found in tissues or within individual biological cells themselves, down to the level of individual biomacromolecules such as DNA. This broadened vision of “hierarchical biomechanics” is the scope of the present subsection.

When dealing with biological tissues, it is well known that the mechanical loading they are subjected to may have a substantial influence on the biological processes responsible for maintaining their integrities, or even for their sizes and shapes [249–253]. In more detail, mechanical loading is usually applied macroscopically, and is somehow transferred to the scale of cellular processes or even below, where corresponding mechanical stimuli influence the occurrence and behavior of biological factors. Hence, multiscale mechanobiology involves, per se, multiscale (poro-)mechanics, on the one hand, and multiscale systems biology, on the other hand. Multiscale (poro-)mechanics has already been dealt with in the previous sections, and it has been successfully applied to bone tissue, in order to estimate mechanical stimuli occurring on the cell scale in response to macroscopically applied mechanical loading; see, e.g., [46,254–258]. These examples are based on or inspired by continuum micromechanics as underlying concept. As an alternative, numerical homogenization approaches similar to those described in Sec. 3.1 have been used for studying bone regeneration guided by tissue engineering scaffolds [85–87,259,260]. A rare, quite mature example of hierarchical mechanobiology, documented in Ref. [255], is illustrated in Fig. 9: The macroscopic (oscillating) loads applied to a bony organ, such as a human femur (see Fig. 9(a)), are first converted into local stresses associated with macroscopic RVEs, one of which is shown in Fig. 9(b). This RVE coincides with that of Level 5 in Fig. 1. The aforementioned stresses are then downscaled to (oscillating) pore pressures in the vascular pores, where they affect the biochemical activities of two types of biological cells: osteoblasts and osteoclasts, see Fig. 9(d). Interestingly, physiological load levels translate into vascular pore pressures of the order of several tens of kilopascals [46]: This is the magnitude to which a large set of biological cells are reacting [261–277]. Further downscaling from the macroscopic RVEs of Fig. 9(b), to the level of lacunar pores (see Fig. 9(c)), delivers pore pressures of similar magnitude, stimulating osteocytes, as seen in Fig. 9(e). In more detail, the increase of lacunar and vascular pore pressures stimulates the proliferation of active osteoblasts. The latter builds up new extravascular bone matrix. At the same time, increasing lacunar and vascular pore pressures suppress the production of RANKL, i.e., of the key factor activating osteoclasts, the cells which dissolve extracellular bone matrix.

It is interesting to note that similarly large pressures occur in the extrafibrillar matrices of soft tissues, as has been quantified by the hierarchical tendon model of Morin et al. [136], also referred to in Sec. 3.5. This holds the promise for wide research opportunities unifying not only hard and soft biomechanics, but also hard and soft tissue mechanobiology.

Biological cells may not only populate pores with extracellular matrices making up tissues (the latter may be categorized as porous solids with largely varying stiffnesses), but they may also govern the behavior of complex fluids such as blood. The latter is made up of red blood cells, white blood cells, platelets, and plasma. From the viewpoint of “fluid micromechanics”, these four constituents would appear as “phases” making up a “fluid material system”. However, the nature of the dynamics of blood (standardly referred to as hemodynamics) rendered particle-based methods [278], rather than semi-analytical solid mechanics methods, as particularly appropriate. Accordingly, biological entities such as cells or platelets are modeled as discrete particles. Moving from the blood itself back to the vessels in which it is contained, particle mechanics has also been beneficial in the context of artherosclerosis, where low-density lipoprotein (LDL) molecules enter into the intima, either through mediation of vesicles taking up LDL or through apoptosis-induced leaky junctions in the endothelium [279,280]. In this context, the so-called discrete particle dynamics (DPD) method turned out as particularly effective when it comes to simulating hydrodynamics of simple and complex liquids on a mesoscopic scale [281]. Details on the mathematical implementation of the DPD approach can be found, e.g., in Ref. [278,282–285], see also Fig. 10. A particularly illustrative example concerns the modeling of thrombosis as the interaction between activated platelets, forming aggregations among themselves and with the adjacent endothelial cells [283], see Fig. 11. Of similar interest is the DPD-based prediction of shear fluid stress-dependent initiation of atherosclerosis by means of binding LDL to the blood vessel walls [279].

We finally turn toward single biological cells and their internal structure: Traditional visco-elastic continuum models have been widely used to quantify cell mechanical properties [286]. Such models are convenient to help parameterize the results of mechanical measurements from techniques such as atomic force microscopy (AFM), micropipette aspiration, microrheology, and parallel plate compression [287]. They offer a straightforward and minimal way to interpret differences between experimental conditions and summarize the results of mechanical measurements. When restricted to the timescale and spatial resolution of a specific measurement, these relatively simple mechanical models are in strikingly good agreement with observed mechanical behavior. For example, the force-indentation response obtained from an AFM measurement typically corresponds well to Hertz's theory for elastic bodies [288]. Under this framework, many model implementations have been introduced to analyze AFM experiments using visco-elastic theory, such as adaptations of Hertz theory for the specific contact geometry [289], models that consider the viscous response of the cell or the distinct mechanical properties of the actin cortex [290–292], and FEM analysis to represent the specific geometry of the cell and the position on the nucleus [293].

However, when generalizing mechanical properties across spatial and temporal scales, visco-elastic or poro-elastic constitutive models invariably break down. Instead, living cells appear like glassy materials that exhibit power law rheology when probed across time scales [294–296]. This property has been attributed to the cells active material properties [297]. At the microscopic scale, the mechanical behavior of a cell is governed by the cytoskeleton, a dynamic reorganizing polymer network that forms mechanically relevant structures at a broad spectrum of characteristic length scales, giving rise to self-similarity and fractals [298]. Furthermore, the cytoskeleton exhibits biphasic rheology, with frequent solid-to-fluid dynamics and even self-organized criticality in active structures such as the lamellipodium [299,300].

To date, the aforementioned properties have complicated homogenization attempts to model the mechanical behavior of intracellular material [301]. Indeed, the absence of a clearly defined microscopic (small) length scale prevents the straightforward formulation of “representative volume elements” in the pursuit of macroscopic constitutive equations. Successful attempts to model the mechanical properties of intracellular material at a cellular scale include coarse-grained representations of acto-myosin polymer networks and continuum models based on hydrodynamic theory of active fluids.

Coarse-grained molecular dynamics models of active polymers aim to provide a computationally minimal representation of the network structure formed by the cytoskeleton [302,303]. Typically, they include, as discrete elements, the filaments, and crosslinkers of acto-myosin or microtubule networks. Contrary to full-atomistic resolved molecular dynamics, they do not include atomic interactions, but instead are based on coarse-grained approximations of thermal forces, internal mechanics, and interaction properties [304]. For example, the mechanical properties of polymers are commonly represented using worm-like chain models, whereas the asymmetry in attachment of myosin motors to actin filaments is modeled using a “ratchet” potential [305,306]. Various computational implementations of these discrete coarse-grained elements exist, such as bead-spring models, rod models, and FE representations of fibers [307]. As a whole, this family of models has been very successful in unraveling the mechanics of specific relevant cytoskeletal structures, such as polymerization in the lamellipodium, the attachment membrane or adhesion complexes, and the formation of stress fibers [308–310]. Moreover, they are particularly useful in studying artificial cytoskeletal structures, such as reconstituted acto-myosin networks, of which the well-known composition and chemical conditions greatly facilitate mechanistic model formulation [311]. Vice versa, such artificial structures provide convenient test beds to calibrate and validate simulations of coarse-grained polymer assemblies. Alternatively, these systems can be experimentally studied at submicroscopic scales, using techniques such as active and passive microrheology [312]. Nonetheless, the applicability of this family of models in the context of a biological cell is limited, not only due to computational limitations, restricting the spatio-temporal scale of the domain, but also due to the inability to exhaustively represent the complex environment in intracellular material.

At the other end, hydrodynamic models of active fluids (or active “gels”) are constitutive continuum models that represent a relevant subset of behaviors of intracellular material [313,314]. In poroviscous theory, pioneered in the context of cells by Dembo and coworkers, the intracellular material is modeled as two coupled interpenetrating fluids, the cytosol and a gel-like cytoskeletal network [315]. By addition of signaling-coupled reaction terms for network assembly and disassembly, these models can simulate various cell dynamical processes. Active nematic theory combines the hydrodynamic equations of an active—typically, contractile—fluid with a polar or apolar director field. In doing so, it is able to capture the intrinsic anisotropy of the underlying network and the spontaneous flow as a result of contractile forces generated by myosin activity [316]. Since they are essentially fluid models, these models are generally unable to describe the elastic mechanical response that would be measured in a typical AFM measurement [317]. However, they provide a very useful tool to help understand the intracellular flows that allow the cell to reorganize itself over longer timescales. For example, simulations using these models have been instrumental in explaining the role of the dynamical organization of actin filaments in furrow formation during development [318,319]. Furthermore, the emergence of spontaneous cell division and cell motility has been observed in models that combine active nematics with poroviscous theory [320].

To address the shortcomings of these limiting conceptual frameworks, various hybrid modeling approaches have been proposed that combine discrete coarse-grained networks with hydrodynamic continuum equations or simplified phenomenological models [321]. These attempts have yielded important insights into specific cell mechanical processes such as bleb initiation, protrusion formation, or cell-matrix interaction [322]. Operating at the cellular scale, these models are used to study mechanisms of active force generation and mechanotransduction. As such, they are suitable for the interpretation of traction force microscopy measurements, in which active forces are quantified by measuring the deformation of a well-characterized substrate. For example, a hybrid model that combines discrete tensional cables in a continuum contractile material was able to relate forces exerted by the cell on the substrate to the generated tension in actin stress fibers [323]. However, these hybrid modeling frameworks are generally limited to specific phenomena, with an associated reduction in spatio-temporal scales. As such, the formulation of a general modeling methodology for describing cell mechanical behavior, based on measurable microscopic parameters and amenable to homogenization for multiscale models, remains a major scientific challenge [317].

A perspective toward overcoming this major challenge recently emerged from a very classical concept in continuum mechanics: the principle of virtual power in its most modern format put forward by Germain [324]: This principle is a natural merger of the laws of equilibrium (or motion in case of inertia forces) with the kinematical characteristics of the investigated mechanical system. Accordingly, it states that the virtual power of internal, external, and inertial forces on arbitrary admissible virtual velocity fields defining the kinematical nature of the considered mechanical system, need always to be zero. Given the current atomistic-to-continuum homogenization challenge, identification of the virtual power performed by equilibrated molecular dynamics-derived atomistic forces in different configurations of the investigated system, on virtual velocity fields characterizing macromolecular kinematics in terms of beam-theory-type continuum notions (stretching, twisting, bending) opened a new avenue to study coupled deformation modes in DNA [325], see Fig. 12: Stretching and twisting of these long double-helical macromolecules are, as a rule, coupled, and the type of coupling depends on the deformational state. This may hold the promise to come up with a new type of structural mechanics which is complex and at the same time stable enough, so as to effectively build up a new, mechanics-driven, addition to the current state of the art in genetics and cell biology. These most probably worthwhile endeavors have been recently coined as “mechanobiome” research [326]. It holds the promise to integrate the vast “computational wisdom” accumulated for DNA and other biomacromolecules [327–329] into (micro)structural mechanics terms, allowing for linking this “condensed wisdom” to the wide field of hierarchical biomechanics. More precisely, any characteristics emerging through homogenization over the atoms making up the biomolecules, and over the latter seen as structural members with appropriately defined mechanical properties, may provide radically new knowledge for at least two challenging open fields in bio-engineering and the life sciences: (i) deeper, and more quantitative, understanding of mechanical stimulation properties of overall biological cells, which could then enter hierarchical models as the one seen in Fig. 9; and (ii) more comprehensive insight into the mechanical competence of extracellular matrices under extreme, i.e., also pathological, conditions, as additional input for hierarchical models as the one shown in Fig. 1.

A set of several scale-specific representative volume elements, nested among each other and each equipped with stress and strain average rules, as well as with scale transition rules for elastic as well as inelastic strains, allows for the prediction of mechanical properties from composition and organization of biological tissues. This has been shown in greatest detail in the case of bone, by combining analytical or computational homogenization schemes with the data obtained from a variety of biochemical, biophysical, and biomechanical experiments. Recent research results suggest a similarly, high potential for the soft tissue regime, and hierarchical modeling approaches seem particularly promising in the fields of mechanobiology (where biological cells are “inserted” into the tissues' pore spaces) and cell/biomacromolecule mechanics.

The authors gratefully acknowledge the financial support granted by the European Commission in the framework of the project Increasing scientific, technological, and innovation capacity of Serbia as a widening country in the domain of multiscale modelling and medical informatics in biomedical engineering (SGABU), Grant Agreement No. 952603.

• European Commission (Funder ID: 10.13039/501100000780).

AFM =

atomic force microscopy

CT =

computed tomography

DPD =

discrete particle dynamics

FE =

finite element

FFT =

fast Fourier transform

ICP-OES =

inductively coupled optical emission spectroscopy

LDL =

low-density lipoprotein

RANKL =

receptor of activator nuclear factor κB ligand

RVE =

representative volume element