1. Introduction

In the present context, soft biomaterials are amorphous solids made of synthetic or natural polymers interacting with or being components of living systems. They may include hydrogels made of synthetic polymers or biological macromolecules (e.g., proteins, polysaccharides), soft tissues or decellularized extracellular matrix (ECM). Structurally, these biomaterials feature cross-linked networks of fibrous constituents, and they exhibit nonlinear and time-dependent behavior in their mechanical response to finite deformation.

That soft biomaterials accumulate enormous attention is attributable to their versatile application in biomedicine—thanks to superior biocompatibility, ease of formation into specific geometries, and broadly tunable physico-chemical properties ^[1,2,3,4,5]. Soft biomaterials often function primarily as scaffolds in artificial living systems; therefore, their mechanical properties directly determine the mechanical characteristics of the whole constructs. Furthermore, recent research has revealed that the mechanical properties of cell-culture scaffolds or substrates are critical physical cues for cellular-fate processes including differentiation of stem cells ^[6,7,8,9]. These discoveries make the study of mechanical properties of scaffold materials and their interactions with cells (i.e., mechanobiology) an incontestably important topic in tissue engineering and stem cell biology.

The mechanical properties of any material, soft biomaterials included, are best understood with the aid of theoretical models. In this review, we seek to integrate some of the major theoretical methods and models used for soft biomaterials and to highlight recent developments in micromechanics, in order to promote a pursuit that is fundamentally important to mechanobiology—the understanding of materials at cellular scale ^[10,11,12]. The wider interest, including such topics as soft biomaterials fabrication, chemistry, rheological experiments, and numerical methods, are served by excellent reviews listed here ^[13,14,15]. In fact, it should be borne in mind that the theoretical methods and models discussed herein are not limited to soft biomaterials, but are generally suited to a broad range of material types.

We begin by reproducing some typical rheological phenomena using the one-dimensional standard linear solid model (SLS model). The SLS model falls into the category of "phenomenological viscoelastic models"—one of the five categories, described at the end of this section, that we use to classify all methods and models discussed in this review. We begin with this model because it is historically significant and incorporates fundamental concepts that illustrate important mechanical behavior encountered in soft biomaterials. Figure 1a is a schematic representation of the model. It can be written in differential form as

$\sigma (t) + \frac{\eta }{{{k_2}}}\dot \sigma (t)\, = \,{k_1}\varepsilon (t) + (\frac{{{k_1}}}{{{k_2}}} + 1)\eta \dot \varepsilon (t)$

(1a)

or in integral form as

$\sigma (t) = {k_1}\varepsilon (t) + {k_2}\int\limits_0^t {{e^{ - (t - \tau )/\lambda }}\dot \varepsilon (\tau )\,d\tau }$

(1b)

where σ and ε are stress and strain, respectively; a diacritic dot indicates the derivative with respect to time t; k₁ and k₂ are the elastic coefficients; η is the viscous coefficient; λ is the relaxation coefficient $\lambda = \eta /{k_2}$.

Figure 1b depicts the stress-strain relationship (SSR) as strain ramping at constant rate. It can be seen that the SSR is influenced by strain rate and that as strain rate approaches zero (i.e., near equilibrium deformation) the SSR becomes the Hookean relationship $\sigma = {k_1}\varepsilon$. Figure 1c shows the stress profile under ramp-and-hold deformation. The stress relaxes exponentially during strain holding. In order to maintain the stress at its maximum in Figure 1c, the strain must continue as in Figure 1d (i.e., the phenomenon of creep).

One common protocol in rheological tests is to impose a small-amplitude oscillation on the sample. In this case, the above SLS equations give rise to a complex modulus:

$G = G' + iG'' = ({k_1} + {k_2}\frac{{{\lambda ^2}{\omega ^2}}}{{1 + {\lambda ^2}{\omega ^2}}}) + i\frac{{{k_2}\lambda \omega }}{{1 + {\lambda ^2}{\omega ^2}}}$

(1c)

where ω is angular frequency; i the imaginary unit. Figure 1e shows the dependence of G' and G'' on ω. The existence of G'', which arises from the phase difference between stress and strain, indicates the existence of viscous dissipation in the material.

These rheological phenomena are often encountered in soft biomaterials. However, the rheological behavior of most soft biomaterials under typical deformational conditions cannot be reproduced with acceptable precision using the SLS model, or a given set of parameters may fit well to one phenomenon and deformational condition, yet fail to represent others. Moreover, there are many phenomena that cannot be reproduced by the SLS model even in a qualitative manner.

Because of the diversity and richness of the real world, various theoretical methods and models have been developed to meet the demands introduced by our utilization of soft biomaterials. We classify these by phenomenological/structural and time-independent/time-dependent categories and mainly confine our discussion to 5 types, as listed in Table 1. This classification is not strict because some methods or models can belong to multiple categories. For example, network models dealing with entropic elasticity ultimately lead to hyperelastic forms, whereas hyperelastic models can also be used to the fibrous constituent of the network models. Nevertheless, this classification helps to separate and elucidate the methods and models proposed in the bulk of the literature.

---

2. Phenomenological Hyperelastic Models

The initiative of hyperelastic models is to address the nonlinear elasticity of materials in the equilibrium state. In contrast to a linear SSR, there are different types of curves in the strain-stress plot (Figure 2): curves that are convex to the horizontal (strain) axis, usually called J-type; curves that are convex to the vertical (stress) axis, herein called ρ-type; and curves that combine both features, called S-type. Interestingly, most of the soft biological materials and tissues have J-or S-type SSR, whereas many synthetic polymeric soft biomaterials have ρ-type SSR. Intuitively, we think this difference may underlie the mechanical mismatch issue ^[16,17,18] at the anastomotic site between host tissue and soft implants, and may even stimulate the scar formation response.

Because we are dealing with equilibrium elasticity, we can assume a strain energy density function W and obtain the stress by differentiating it with respect to the deformation tensor. Details of the hyperelasticity and strain energy function may be referred to ^[19,20,21]. Here we list some commonly used models and discuss their roles in the recent literature.

Let F be the deformation gradient tensor, with principal stretches denoted by γ₁, γ₂, and γ₃; the principal invariants of its corresponding left Cauchy-Green tensor are

$\left. \begin{array}{l} {I_1} = \gamma _1^2 + \gamma _2^2 + \gamma _3^2\\ {I_2} = \gamma _1^2\gamma _2^2 + \gamma _2^2\gamma _3^2 + \gamma _1^2\gamma _3^2\quad \\ {I_3} = \gamma _1^2\gamma _2^2\gamma _3^2 \end{array} \right\}$

(2)

The strain energy density functions of neo-Hookean, Mooney-Rivlin, Ogden, and exponential models for incompressible materials are, respectively:

$\left. \begin{array}{l} {W_{NH}} = c({I_1} - 3)\\ \\ {W_{MR}} = {c_1}({I_1} - 3) + {c_2}({I_2} - 3)\\ \\ {W_O} = \sum\limits_{p = 1}^N {{\mu _p}(\gamma _1^{{\alpha _p}} + } \gamma _2^{{\alpha _p}} + \gamma _3^{{\alpha _p}})/{\alpha _p}\\ \\ {W_{\exp }} = {b_1}(\exp ({b_2}({I_1} - 3)) - 1) \end{array} \right\}$

(3)

where c, c₁, c₂, b₁, b₂, μ_p, and α_p are coefficients. Then, the 1st Piola-Kirchhoff stress tensor can be obtained from ${\rm{T}} = \frac{{\partial W}}{{\partial \, {\rm{F}}}}$.

Some modified models and other possibly usable models have been proposed ^[22].

The application of hyperelastic models to biological tissues has a long history ^[20,23]. These models can reproduce all of the aforementioned types of SSR (J, S, and ρ) for soft tissues. Recently, Madireddy et al. ^[24] employed Mooney-Rivlin, Ogden, and exponential models to investigate the nature of agarose gel, bovine liver tissue, and porcine brain tissue, and showed that the two-term and three-term Ogden models provided better fits than did the Mooney-Rivlin and exponential models. In particular, they developed a Bayesian calibration framework based on nested Monte Carlo sampling for parameter estimation, rather than the traditional least squares method, to address the model identification issue for choosing between hyperelastic models.

Biological tissues are generally anisotropic, and this is mainly attributed to constituent/fiber orientation and fibrous crimp. Therefore, researchers have for some decades incorporated these structural elements when applying hyperelastic models to biological tissues ^{[25,26,27,28]}. In dealing with the cyclic hysteresis and time-dependent phenomena it has been realized that structural alteration must be taken into account, which introduces the damage function into the modeling ^[29,30,31]. These modifications, which are based on structural considerations, feature the application of phenomenological hyperelastic models to biological tissues and will continue to develop.

Recently, Horgan and Murphy ^[32] conducted a theoretical investigation of shearing response in anisotropic soft biological tissues that incorporated these modifications. Upon reproducing the J-type shearing response in biological tissues, they pointed out ambiguity in the determination of normal stresses and suggested an accurate shear test protocol for the constitutive properties of soft tissues. Amar et al. ^[33] tested the fibrous capsular tissue surrounding breast implants by uniaxial tension. They found that it was difficult to fit the data over the entire range of tested strain with a single hyperelastic model, so they used the Mooney-Rivlin model at low strain and the Valandis-Landel model at large strain. It was revealed that the models were too sensitive to fibrous orientation to be able to identify the parameters. These studies raise a critical issue in the choice of a hyperelastic model for anisotropic soft biomaterials: there is no universally applicable hyperelastic model. Verification and parameter identification for model choice usually require large deformation and multi-dimensional tests, which are difficult to conduct with precision on soft tissue samples.

In the area of micromechanical testing on soft biomaterials, Lin and Horkay ^[10] have surveyed progress in the application of atomic force microscopy (AFM) and nanoindentation technology to characterize the local mechanical properties of polymer gels and biological tissues, in consideration of various hyperelastic models. They emphasized the need for multiple models to span the diverse behavior of soft biomaterials. The authors of this review question the applicable regime of the hyperelastic model in micromechanically dealing with the manifest structure of soft biomaterials. Hyperelastic modeling is based on continuum mechanics and aims to reveal the macro response of materials. At how small a scale the theory still holds is an issue deserving of careful investigation.

A further application of hyperelastic models in micromechanical characteristics of soft biomaterials is the description of the mechanical properties of a single composition element within structural models, such as in biphasic or network models; this will be returned to in sections 4 and 5.

---

3. Phenomenological Viscoelastic Models

As with models driven by phenomenon and experimental data fitting, this family of models does not attempt to explain the causal mechanism of time-dependent mechanical response; therefore, the meaning of viscoelasticity in these models is different from that in the structural models (described in the next section), where time-dependent behavior is distinguished by "viscoelasticity" intrinsic to the solid structure and fluid permeation within the structure called "poroelasticity". The viscoelastic property of phenomenological models (in this section) combines these phenomena.

The SLS model, mentioned in the Introduction, is representative of this category. Although the validity of such linear viscoelastic models is limited, this type of model is still a good choice owing to its linearity and the simple form for synthetic polymers under small deformation. Saxena et al. ^[34] tested agarose, poly-2-hydroxyethylmethacrylate gels, and rat spinal cord tissue using their microindentation system and found that the stress relaxation behavior in these materials could be regressed via the SLS model with good agreement. In ^[35], White et al. comprehensively described the traditional and nano-indentation technology; they employed Kelvin-Voigt viscoelastic model to represent the viscoelastic characteristics of the instruments and material samples. By testing four different synthetic polymers they found good agreement between nano-indenter data and data produced by a traditional analyzer, in materials with weak viscous features.

When applying this kind of classical method to soft biomaterials, issues encountered during experimental parameter identification become critical. Due to the softness and labile nature of the materials, experiments are sensitive to the sample setting and the initial conditions of the test. In this regard, Tirella et al. ^[36] proposed a method (the "epsilon dot method") that eliminates the influence of pre-strain/pre-stress during sample installation when deriving the viscoelastic parameters in SLS or 5-parameter models characterizing the properties of hydrogels of polydimethylsiloxane (PDMS) and gelatin, and fresh porcine livers. The swiftness and precision of the method for parameter identification in stress relaxation and creep experiments were proved.

In dealing with soft biomaterials composed of natural components such as collagen, fibrin, actin, or ECM extracts, however, modeling these materials becomes complicated due to the strong nonlinear properties ^[37,38]. A straightforward approach to this is to nonlinearize the element in linear models, which includes nonlinearzing the elastic element by using hyperelastic models ^[39,40] and making the relaxation coefficient a function of strain ^[40,41]. In this respect, Fung's quasi-linear viscoelastic (QLV) model is well-known ^[21]; it can be written as follows:

$\sigma (t) = Q[\varepsilon (t)] + \int\limits_0^t {Q[\varepsilon (t - \tau )]\frac{{\partial G(\tau )}}{{\partial \tau }}\,d\tau }$

(4)

where the function Q[ε(t)] is called the elastic response and G(τ) the reduced relaxation function. The schematic representation of this model is shown in Figure 3a. This model asserts that the stress at any time t equals the material elastic response decreased by some quantity that depends on the history of this same elastic response, weighted by the time differential of the reduced relaxation function. This model finds good application to biological tissues including vessels, tendons, ligaments, cartilage, and muscle.

Recently, we proposed a non-linear Kelvin-type model with two relaxation coefficients ^[42]:

$\sigma (t)\, = \,{k_1}{\varepsilon ^\alpha }(t) + {k_2}\int\limits_0^t {[\xi {\kern 1pt} {e^{ - (t - \tau )/{\lambda _1}}} + (1 - \xi ){\kern 1pt} {e^{ - (t - \tau )/{\lambda _2}}}]} \frac{{\partial {\varepsilon ^\beta }(\tau )}}{{\partial \tau }}d\tau$

(5)

The idea in this model is that external work is split into two fractions (Figure 3b): the first is reversible energy storage in a power-law spring, with coefficient k₁, and the second is a power-law spring with coefficient k₂, in which energy is dissipated by two exponential processes with relaxation coefficients λ₁ and λ₂ (λ₁ < λ₂). This model can fit data obtained by uniaxial tensile and compression experiments on decellularized ECM, ECM extract hydrogels, and biological tissues ^[42,43,44]. It can also reproduce the different patterns of the dependence of the elastic reversible fraction on the holding strain ^[42].

The relaxation process of most soft biomaterials possesses a fast phase and a slow phase. This can be reproduced using a pair of exponential functions with small and large relaxation coefficients, as in Eq. (5). The problem here is that the stress will eventually decay to a constant value determined by the larger coefficient λ₂. For a mechanical process that lasts no longer than a few multiples of the coefficient λ₂, a good fit can always be achieved. However, in some soft biomaterials, relaxation can continue for days without settling at a constant value. One method for tackling this problem is to parallelize many more viscoelastic models, each with different relaxation coefficients; however, this method is obviously less practical in parameter identification.

A power-law relaxation might provide a better treatment to this issue. The introduction of a fractional derivative into the viscoelastic constitutive law was pioneered by Scott-Blair ^[45]. Later, Schiessel et al. used a fractional element in traditional linear viscoelastic models and obtained closed analytical solutions showing power-law relaxation with respect to time ^[46]. Recently, Jaishankar and McKinley ^[47] developed a fractional Maxwell model (Figure 3c) that was highly effective in reproducing power-law relaxation, complex modulus behavior, and creep compliance on butyl rubber, acacia gum, and bovine serum albumin interfaces. de Sousa et al. employed a fractional SLS model to investigate the viscoelasticity of polyacrylamide gels under AFM ^[48]. They found that with increasing bisacrylamide concentration the viscoelasticity of the gels became more and more significant. The potential of this kind of model for characterizing soft biomaterials deserves considerably more attention.

---

4. Structural Biphasic Models

The progression of modeling approaches from phenomenological to structural is effectively an evolution in understanding from "how it is" to "what it is".

Because the primary structural feature of soft biomaterials is their high water content, biphasic models are the natural choice for mechanical characterization. Biphasic theory was created by Biot to treat the properties of soil ^[49] and extended by Bowen ^[50], Green and Naghdi ^[51], and Mow ^[52] to biological tissues. The core concept of the biphasic theory is as follows. If a mixture material is supposed to be constituted by a solid skeleton and an interstitial fluid, one can obtain the continuity equations and equations of motion, respectively, for each constituent and for the overall mixture. The coupling of these two constituents appears in the equations of motion as a term known as "diffusive force". These governing equations must be enclosed by constitutive equations for each constituent, with an additional expression for the diffusive force. In the framework from Bowen's to Mow's, the solid skeleton can be elastic or linear viscoelastic, the fluid is viscous, and the expression for diffusive force is obtained by thermodynamic consideration of the viscous dissipative flow. In the original framework of Biot's, linear deformation approximation is adopted, the solid is elastic, the fluid is inviscid, and the diffusive force becomes an incremental hydrostatic pressure related to deformation through Darcy's diffusion law. From this core concept, there have been a number of theoretical modifications corresponding to particular concrete applications ^[53].

The advantage of the application of biphasic models to soft biomaterials is that they are particularly apt for modeling the behavior under compression, while the mechanical properties are directly related to the physical and chemical parameters of the materials.

Based on Mow's biphasic framework and a fluid transport model suitable for biological soft tissue, Gao and Gu ^[54] proposed a new model for soft tissues and hydrogels and investigated the effects of hydration on mechanical properties including Poisson's ratio. The predictions of the aggregate modulus and shear modulus for human and animal cartilage and for poly (DMAEMA-co-AAm) hydrogels and agarose hydrogels by the model compare well with those from experimental results. Most recently, Bonilla et al. ^[55] modeled the irrecoverable compression of cellulose/xyloglucan composites by introducing an aggregation force and a critical yield pressure into the constitutive biphasic formulation for transversely isotropic tissues.

Because the structural feature is captured in biphasic models, biphasic constitutive theory is able to explore the origin of viscoelasticity and extends to the discrimination between viscoelasticity in the solid skeleton and poroelasticity due to fluid permeable movement. By employing Fung's QLV model for the solid phase of the agarose gels in a biphasic model, Olberding and Suh numerically simulated creep and relaxation processes recorded by indentation experiments ^[56]. They found that agarose gels possess intrinsic viscoelasticity in the solid skeleton—without consideration of which the creep and relaxation processes could not be reproduced. In the case of collagen gels, however, Castro et al. ^[57] showed by numerically simulating the compression-relaxation process within the biphasic framework incorporating a neo-Hookean solid phase that the rheological behavior of gels under compression was mainly determined by poroelasticity. Interestingly, they reproduced fast stress decay, which is usually regarded as a consequence of structural viscoelasticity.

In general, the time-dependent behavior of soft biomaterials is attributable to the co-existence of structural viscoelasticity and permeable poroelasticity ^[58]. In stress relaxation processes, structural viscoelasticity manifests as fast stress decay with a large fraction in the decreased stress, with cross-link concentration dependence and sample size independence; conversely, permeable poroelasticity manifests as slow stress decay with a small fraction in the decreased stress, and is typically dependent on sample size. Biphasic theory has also been applied to analyze the micro-indentation experiments. Studies showed that poroelasticity in synthetic hydrogels resulted in larger modulus in micro-indentation experiments than that in macroscale homogeneous ones and confirmed that the poroelastic time-dependent behavior was length scale related ^[59,60]. The theory also showed that polyacrylamide gel—usually regarded as a pure-elastic, time-independent material—actually exhibited stress relaxation in nano-indentation testing but not in micro-indentation testing ^[61]. Wang et al. proposed a method to separate viscoelasticity and poroelasticity by testing samples with sizes larger or smaller than corresponding material lengths ^[62].

---

5. Structural Network Models

The motivation behind structural network models is to account for the mechanical properties of rubber materials ^{[63,64,65,66]}. The statistical study of the mechanical response of single synthetic polymer chain lays the cornerstone of this family of models. These polymer chains act as flexible entropic elastomers whose force response is determined by temperature and the end-to-end distance of the chain. Depending on the polymeric species, the statistical distribution of the end-to-end distances can be Gaussian or non-Gaussian. When determining the bulk mechanical properties by integrating individual force response, the network models come into play. Therefore, network models in the continuum regime actually include two parts: the force response of the individual chains and the network itself, which connects the end-to-end distance to the macroscopic deformation.

The basic framework of a network model may be expressed as follows. The Helmholtz free energy of single chain is

$\Psi = E - \theta \,S$

(6)

where E is enthalpy energy; θ is absolute temperature; and S is entropy. For a flexible chain, E can be neglected; then, by single-chain statistics, the entropy $S(e) = {k_B}\ln P(e)$ can be found, where e is the end-to-end distance, P(e) is its distribution function, and k_B is Boltzmann's constant.

To obtain the macroscopic mechanical constitutive relationship, we need to differentiate the Helmholtz free energy of the whole bulk with respect to strain F

${\rm{T}} = - \theta \,\frac{\partial }{{\partial {\rm{F}}}}\int {\int {p(e)S(e)} \;de\,dV}$

(7)

where p(e) is the distribution of e in dV. The traditional network model assumes there are μ (number) chains in dV, each chain having the same end-to-end distance ${\bar e}$. Here, the goal of the network model is to connect ${\bar e}$ with F:

$\bar e \Leftrightarrow {\bf{F}}$

(8)

For example, the well-known 8-chain model ^[67] sets out that

$\int {P(e)S(e)\;de = S(\bar e):\,\quad \quad \bar e \Leftrightarrow \frac{a}{2}} {(\gamma _1^2 + \gamma _2^2 + \gamma _3^2)^{\frac{1}{2}}}$

(9)

where a is the dimension of dV; γ₁, γ₂, and γ₃ are the principal stretches (Figure 4a).

In addition to the above method, which integrates a geometric network model into a continuum scheme, another class of discrete network models, which are artificially generated via a particular algorithm or pattern (Figure 4b), also plays a critical role in structural network theory. The discrete network models can be further categorized as lattice-based models ^[68] and off-lattice models ^[69]. These models are particularly suited to investigating the mechanical properties of biological hydrogels at low constituent concentrations. The advantage of this class of models is that it is easy to reproduce many features of soft biomaterials, including non-affinity, heterogeneity, and anisotropy.

The network model found ready application in soft biomaterials because most of these are of a hydrated fibrous structure. However, we would like to emphasize that unlike rubber materials many soft biomaterials are composed of semiflexible chains; these are distinct from flexible chains in that: (1) the statistical treatment of S(e) is much more complex, depending on the contour length and persistence length of the chain; and (2) the enthalpic energy E in $\Psi$ cannot be neglected but rather may dominate $\Psi$, depending on the molecular species or deformation. Broedersz and MacKintosh have thoroughly reviewed progress in semiflexible polymers and network models ^[70]. The authors would like to highlight several excellent works on the statistics of semiflexible chains ^[71,72], chains possessing both enthalpic and entropic elasticity ^[73,74,75], and networks of semiflexible chains ^[76,77,78].

Recent studies continue to elucidate the important roles of network characteristics in determining the mechanical properties. Experimental data are further enriching our understanding of network characteristics including non-affine deformation ^[79], heterogeneity ^[60,80], anisotropy ^[33], and structural hierarchy ^[81]. These features challenge traditional networks based on homogeneous and isotropic assumptions. One way to resolve this disparity is to numerically generate artificial networks ^{[68,69,82,83]} and embed these features.

For theoretical methods, Cioroianu et al. have modified the classic 8-chain model to assess the role of certain network distortions, such as disorder, pre-stress, and non-affinity ^[84]. By off-positioning the central node of the 8-chain network, they found that even when each individual chain behaves linearly, the whole polymer model exhibits nonlinear behavior. Actually, it has been demonstrated through discrete network models that nonaffinity of the filamentous network can lead to strain-stiffening of soft biohydrogels ^[85,86].

Recently we proposed a new network framework aimed at coping with features of real fibrous materials ^[87]. The basic concept can be roughly explained as follows.

Considering that a semiflexible chain may enter the enthalpic regime, we keep the enthalpic term in the expression of Helmholtz free energy and suppose it can continue to be expressed as a function of e, thus:

$\begin{array}{l} {\rm{T}} = \frac{\partial }{{\partial {\kern 1pt} {\rm{F}}}}\int {\int {p(e)\Psi (e)} } \;de\,dV\\ \quad = \int {\int {\,(\frac{{\partial p(e)}}{{\partial {\rm{F}}}}} } \Psi (e) + p(e)\frac{{\partial {\kern 1pt} \Psi (e)}}{{\partial {\rm{F}}}})\;de\,dV \end{array}$

(10)

The first term in the integrand represents the distribution change due to deformation. It can be classified into network deformation and crosslink breakage. We assume that the contribution of network deformation can be neglected comparing to that of crosslink breakage. Thus, this term addresses the network alteration characteristics, which will be explained further in the next section. The second term can be further expressed as:

$p(e)\frac{{\partial {\kern 1pt} {\kern 1pt} \Psi (e)}}{{\partial {\rm{F}}}} = p(e)\frac{\partial }{{\partial {\rm{F}}}}(\Psi ({e_{iso}}) + \Psi ({e_{ori}})) = p(e)\frac{{\partial {\kern 1pt} \Psi }}{{\partial e}}(\frac{{\partial {e_{iso}}}}{{\partial {\rm{F}}}} + \frac{{\partial {e_{ori}}}}{{\partial {\rm{F}}}})$

(11)

where $\frac{{\partial \Psi }}{{\partial e}}$ is the force response of a single chain; e_iso and e_ori are those e with random orientation and those with a specific orientation, respectively (Figure 4c). e_iso is an isotropic function of F, thus have to be under the argument of the invariant of F; e_ori can be obtain by affine condition. The details may be referred to ^[87]. Obviously, for 8-chain model ${e_{iso}} = \frac{a}{2}{(tr{{\rm{F}}^2})^{\frac{1}{2}}}$ and without e_ori; whereas, for 3-chain model ${e_{iso}} = \frac{a}{3}(tr{\rm{F}})$.

Therefore, this theoretical framework allows us to investigate the effects of end-to-end distance distribution and the effects of fibrous orientation (i.e., anisotropy) on mechanical properties. Moreover, our theory distinguishes itself in that it constitutes a complete theoretical framework for structural modeling, being composed of the force response of a single constituent (including entropic and enthalpic elasticity), the above described network model, and the network alteration theory explained in next section. This theoretical framework can reproduce and predict a wide range of nonlinear mechanical behavior in soft biomaterials, including stress relaxation, hysteresis under serial cyclic loading, strain-stiffening, and the so-called "negative normal stress" phenomenon. Applied to compacted collagen gels, the theory demonstrates that collagen fibrils behave as enthalpic elasticas with linear elasticity within the gels, and that the macroscale nonlinearity of the gels originates from the curved fibrillar network.

---

6. Structural Alteration Theory

Soft materials, including soft biomaterials, often experience large cyclic deformation, consequently exhibiting hysteresis and the Mullin effect—instantaneous and irreversible softening of the stress-strain curve under continuous cyclic loading-unloading with incremental maximum strain in each cycle. It is difficult to reproduce such phenomena theoretically without taking into account structural alteration.

Structural alteration theory was originated by AV Tobolsky who proposed that a portion of the initial network of polymeric materials is broken and reformed into a new network ^[88]. The concept is illustrated in Figure 5. The unaltered equilibrium portion can be modeled by entropic elastic models or phenomenological hyperelastic models. The question is how to deal theoretically with the altered portion and the network evolution to relate it to the prior deformation. Tobolsky's method ^[89], based on the crosslinking chemistry, is to introduce an exponential function of deformation time t for chain breakage,

$s = {s_0}{e^{ - kt}}$

(12)

and a chain-growth function for network reformation in the time interval $dt'$,

$ds = skdt'$

(13)

where s represents the instantaneous number of chains; and k is consequently the relaxation coefficient. The important point in this theory is that broken chains instantaneously enter a stress-free state, whereas the reformed chains develop a stress that is determined at the moment of reforming nucleation.

In the continuum mechanics scheme, the structural evolution is evaluated by a damage function controlled by the prior maximum strain ^{[29,30,31,90,91,92]}; for instance, it can be expressed as

$g({\varepsilon _{\max }}) = \beta + (1 - \beta )\frac{{1 - {e^{ - \,\frac{{{\varepsilon _{\max }}}}{\alpha }}}}}{{\frac{{{\varepsilon _{\max }}}}{\alpha }}}$

(14)

where ${\varepsilon _{\max }}$ is a kind of maximum deformation; α and β are material constants. The damage function is essentially phenomenological and empirical, and deformation time t does not explicitly appear.

Inspired by the reptation of macromolecules (entangled or with free ends), Bergstrom and Boyce ^[93,94] proposed a formulation to relax the altered structure leading to time-dependent behavior, and applied this to biological tissues.

Development of the alteration formulations, as evidenced by many detailed structural concerns corresponding to specific materials, is continuing. In consideration of the cross-link as being trapped in a potential energy well, Meng et al. ^[95] proposed expressions for cross-link breakage and reformation based on Kramers' kinetics. Li et al. ^[96,97] tackled the elastomer network alteration issue under large deformation by using a modified reptation tube model and incorporating non-affine deformation. To model the superior mechanical properties of nanoparticle crosslinked hydrogels, Wang and Gao ^[98] developed a network model composed of polymer chains with inhomogeneous length in the 8-chain network crosslinked by nanoparticles. Chains are assumed to detach from the nanoparticles during loading, given that the chain forces are greater than the bonding strength between the particle and the polymer chain; chains may re-attach during unloading. The rates of chain detachment and re-attachment obey exponential formulas, under arguments based on the maximum of the first invariant of the deformation.

To explain the recently discovered strain-enhanced stress relaxation in weakly cross-linked biological hydrogels, Nam et al. ^[99] employed Bell's dynamics of force-dependent binding probability ^[100] to predict crosslink breakage in an image-based network model; therein, an interesting probability rule for determination of the breakage via binding probability was applied. In contrast to the dominance of crosslink breakage under tensile or shearing deformation, crosslink formation is highly significant during compression. In ^[101], Vos et al. simulated the compressive process of fibrin gels by using a 2D lattice-based network model and presented, experimentally and theoretically, deep insight into the compressive stiffening mechanism of filamentous soft biomaterials. This study is particularly relevant to the fabrication of tissue-engineered constructs by means of plastic compression of biohydrogels ^[102,103].

In ^[87], we tackled network alteration by reforming the first term in the integrand in Eq. (10) as an energy dissipation process due to cross-link breakage. Departing from traditional structural alteration theory, which instantaneously sets the broken portion of the network to the stress-free state, we modeled the process of post breakage by a mass-spring-dashpot system and showed that force relaxation after breakage is intrinsically an exponential decay process with two relaxation coefficients. This method is significant in understanding the origin of viscoelasticity in soft materials, and precisely reproduces the hysteretic and Mullin-like phenomena observed in cell-compacted collagen gels.

---

7. Examples of Application to Cell Biology

With the progress of mechanobiology and its potential in tissue engineering and stem cell engineering, the mechanical interaction between the cell and its surrounding biomaterial environment and the consequent effects of this interaction on cell biology are accumulating enormous research interest. We introduce some examples of the application of theoretical methods and models in this field.

Tissue-equivalents fabricated by embedding fibroblasts into collagen gel have played an important role in understanding cell-fibrous biomaterial interaction. The contraction of collagen gel due to this interaction was first described in detail by Bell et al. ^[104]. Since then, it has been regarded as the standard tissue model in studies of the wound healing process and cellular fate processes, and in applications in regenerative medicine ^[105]. It is particularly attractive for mechanical explanations of the phenomenon of gel contraction. Barocas and Tranquillo ^[106] were the first to provide sophisticated mechanical modeling of this issue. They regarded collagen gel as a continuum biphasic material wherein the cells are a component of the network phase, which is phenomenologically modeled as a Maxwell-like fluid. Here, collagen gel contraction involved cell-generated stress, cell concentration, and cell orientation. Zahalak et al. ^[107] developed a constitutive model for capturing the mechanical properties of the construct. They represented the collagen matrix as a QLV continuum and separated cell function using anisotropy tensors to consider the effects of cell orientation. Within this combination of continuum and statistical mechanics, cell traction force and cellular passive mechanical properties could be detected.

We studied the phenomenon from the perspective of collagen fibril mechanics ^[108]. We classify the collagen fibrils into three types—bent, stretched, and adherent—and deduce the respective equations governing the mechanical behavior of each type. Via careful verification of a structural elementary model based on this classification, we paint a clear physical picture of the contraction process, quantitatively elucidate the panorama of the micro mechanical niche, and reveal an intrinsic biphasic relationship between cellular traction force and matrix elasticity.

Recently, in addition to the aforementioned averaging approach, modeling of the whole construct and cell behavior through single cell analyses is becoming prominent. This trend is promoted by progress in experimental technology ^[109,110] that allows inhomogeneity, anisotropy, and asymmetry around individual cells to be measured.

Zeng and Li ^[111] developed a multiscale soft matter model to help understand the mechanotransduction mechanism. Here, the cell was modeled as a hyperelastic medium wrapped by a nematic liquid crystal (i.e., the cell membrane model) while the substrate was hyperelastic. Contact was modeled as an attractive potential force counteracted by a repulsive force calculated via continuum contact mechanics. Using this model, they investigated the effect of substrate elasticity on cell adhesion. To investigate the propagation of mechanical signals through hydrogels, Aghvami et al. ^[112] modeled the substrate as a hyperelastic continuum, a fibrous growth network, and a fibrous Delaunay network, respectively. Cell traction force was uniformly distributed on a local area of the substrate grid. It was found that cell traction force penetrated farther into substrates with low fiber connectivity.

The cell is itself a soft biomaterial. Representing a cell using the methods and models reviewed here is a natural option. For example, cytoplasm can be modeled as a viscoelastic continuum ^[113,114] or as a biphasic material with poroelasticity ^[115,116]. More literature on this topic can be found in ^[115,117].

---

8. Challenges and Perspectives

Table 2 summarizes the methods and models introduced in this review with their pros and cons noted, along with landmark references and recent progress.

Hyperelastic models encounter a dilemma upon application to soft biomaterials. On one hand, hyperelastic models were invented to deal with nonlinear materials under finite deformation in the equilibrium state, which is usually understood as testing at an extremely slow deformation rate. On the other hand, structural alteration of soft biomaterials under finite deformation is inevitable no matter how small the deformation rate ^[29,30,31]. Thus, this dilemma limits their application to soft biomaterials under finite deformation. Another issue is that the continuum approach to hyperelastic models demands that the materials be homogeneous and isotropic, which is obviously not true of most soft biological materials ^{[25,26,27,28]}. These cons compromise the ability of hyperelastic models to deal with nonlinearity of soft biomaterials under finite deformation, and limit their suitability to the realm of homogeneous synthetic flexible polymers or to the description of the mechanical properties of the solid constituents in a structural model.

Phenomenological viscoelastic models have a history as long as that of continuum mechanics, but the challenges of applying these sophisticated models to soft biomaterials are more experimental than theoretical ^[32,33,118]. Due to the anisotropic, heterogeneous, soft, and labile nature of most soft biomaterials, the results are extremely sensitive to the initial and boundary conditions in the experiments. This makes parameter identification through experimentation a tough issue even with careful consideration of the sample setting ^[36,24]. One solution to this issue is to adopt guidelines for standardizing in vitro test conditions—including equipment, temperature, strain rate, sample initial state, and deformational process—to make them relevant to in vivo conditions, or to adopt the conditions suggested by Mattei and Ahluwalia in a review of the mechanical testing of liver tissue ^[118].

This difficulty in application of phenomenological models to soft biomaterials, together with the need to answer intriguing questions raised in mechanobiology, catalyze the trend of structural methods and models in the field. To devise an unambiguous physically based model is, however, still fraught with challenge. We examine these challenges over the development of network models of soft biomaterials.

First, we look into the force response of a single constituent chain in network models.

The commonly used models in the literature for single-chain soft biomaterials are: 1) Langevin statistic model ^[67,92,93], which strictly speaking is only suitable for entropic flexible chains; 2) the well-known asymptotic formula obtained by Bustamante et al. ^[119] through stretching of DNA ^[71], which is suitable for entropic semiflexible chains with contour length much longer than persistence length; 3) an approximate expression based on statistical mechanics taking into account the stretching and bending energy of the chain ^[69], which is suitable for short semiflexible chains deformed in both entropic and enthalpic regions; and 4) our recently deduced empirical formula for semiflexible chains under entropic-enthalpic deformation ^[87], based on the work of ^[71] and ^[75].

Obviously, the first and second single-chain models are not compatible with most soft biomaterials made of semiflexible chains with short network meshes. However, a number of studies that employed these models obtained good agreement with experimental data on bulk mechanical properties; the reason will be explained below. The third and fourth models are substantially universal in their handling of diverse deformations, making them suitable for the materials of interest. However, owing to the lack of direct measurements of the material parameters for these models, these parameters can only be inferred by fitting experimental data to bulk mechanical properties. Due to the variety of the network models, choosing a correct model for individual chains is essential for a physically based method. Therefore, there is a need to establish a complete database of mechanical tests of the most commonly used macromolecules and fibers in soft biomaterials, including experimental data on persistence length, force response at different contour lengths and end-to-end distances, and the transition of stretching-bending for short semiflexible chains. This may be a difficult task, but it is entirely feasible given current experimental and theoretical progress in the micromechanics of single macromolecules and fibers ^{[120,121,122]}.

Secondly, we look into the network models.

The physical manifestation of the elementary constituent of any network model is a fiber segment sectioned by two neighboring crosslink points. To describe this elementary constituent, we need two sets of physical parameters: one for the status of the segment, including end-to-end distance (crosslink density), contour length, and orientation; and the other for the status of the crosslink, including crosslink strength and torsion (or absence of torsion). In consideration of the deformation, one more important condition is the affine or non-affine deformation of the constituent to the macro-deformation. Therefore, ideally, a complete network model should consist of all the above parameters based on experimental measurement together with a correct single-chain response also experimentally determined. However, in practice, due to lack of the direct experimental data in the micro scale, most of the parameters are regarded as variables to be determined within a reasonable physical range by fitting the bulk experimental results. Our opinion is that good agreement with the bulk experimental results does not ensure that a model is physically correct.

Consider collagen hydrogel, for example. Recent studies, including ours ^[87,123,124] have shown that collagen fibrils in collagen hydrogels behave as athermal elastica, and that the nonlinear strain-stiffening behavior of the gels is due either to the non-affine deformation in gels with low density collagen or to the crimped fibrils in fibroblast-compacted gels with high consequent collagen density. However, a number of studies using the Langevin statistical model with the affine deformation assumption to describe collagen fiber/fibril in the network model were able to obtain good agreement with the bulk experimental data. We explain this outcome as follows.

Suppose that the constituent element in a network model is not a single segment sectioned by two neighboring crosslink points on a fibril, but rather is composed of multiple such segments as shown in Figure 6. These multiple segments would render the constituent element longer and the deformation of the element more affine with regard to the macro-deformation than for a single segment. Moreover, the multiple crosslink points on this element fibril could behave approximately like random thermodynamic fluctuations. As a result, the entropic Langevin statistical model under the affine deformation condition could accord with experimental data. The parameter for the constituent element length would undoubtedly fall into the physically reasonable range whether it represented one or several segments. Therefore, the current criteria for evaluating a network model are insufficient to ensure that it is physically realistic.

The above analysis also triggers our curiosity: would it not be feasible to establish the kind of network model shown in Figure 6 so long as a sound statistical mechanical description is established for the constituent element composed of multiple segments? This could offer advantages for the appropriate treatment of non-affine deformations, which is currently a difficult problem for the continuum network models.

Compared to continuum network models, discrete network models are superior in dealing with network disorder, such as anisotropy and heterogeneity. In particular, discrete network theory introduced network connectivity ^[68,125] and bending torsion at crosslinks ^[68] as fundamental parameters of the network models to theoretically investigate the mechanical properties of semiflexible biological hydrogels with low content density. Network models with these features are able to interpret some critical phenomena of low density biological hydrogels, such as structural stability, strain/stress-dependent properties, and scaling in the relationships between properties and deformation ^[123,124]. The roles of these network features in determining the properties of high density soft biomaterials are worthy of further investigation. On the other hand, the role played by water in low density hydrogels cannot be ignored and should be investigated ^[57]. Actually, according to our experimental experience, low density collagen hydrogels are unstable, and slowly collapse under a slack state.

In view of the insufficiency of micromechanical information on the constituents of structural models, adding more macro verifiable properties would be a realistic way to establish a sound physically based structural model. We recommend the following two properties:

1) Temperature response

Temperature affects a spectrum of physical parameters in soft biomaterials, making temperature-response experiments complex and the testing of any individual parameters difficult. However, other macroscale phenomena occurring with temperature change, such as phase transition and the pattern change in mechanical property-deformation plots can be predicted by means of refined structural models ^[126] and can be used in turn to validate the models themselves by macroscale experiment. Moreover, our study ^[87] showed that temperature response over different strain ranges was particularly sensitive for certain physical characteristics and could be used to verify the model.

2) Poisson ratio

The Poisson effect is an intriguing issue in soft biomaterials ^{[54,79,127,128,129]}. The Poisson ratio of soft biomaterials is much higher (typically 0.8~4.0) than that of traditional engineering materials (typically < 0.5). However, there is as yet no conceptual theoretical model addressing this issue. We expect refined structural models to reveal the underlying process, and the elucidation of this process to provide a test of the correctness and power of the models.

Compared with the structural network models, structural alteration theory is much more phenomenological and empirical, despite the two theories often being used together in modeling soft biomaterials. The weakness in the detailed physics in structural alteration theory also compromises the physical soundness of those combined theories. The pressing issue here is to establish an in-situ wet experimental system to observe the alteration process with deformation; this kind of observation is technically feasible ^[79,99,130]. Also, the instantaneous stress-free treatment of the aftermath of crosslink breakage in these theories is obviously oversimplified and needs to be refined.

Mechanobiology has revealed that the elasticity ^[6,7] and, most recently, the viscoelasticity ^[8,9] of the cell substrates are able to control cellular-fate processes including differentiation. This paves a promising path to harness physical cues for tissue therapy and regeneration. However, there is a "language gap" hindering the elucidation of the underlying mechanism. This is because concepts such as elasticity or viscoelasticity belong to macroscopic continuum mechanics, whereas the cell does not dwell in the realm of continuum mechanics but can only sense microscopic, discrete, local effects. This is the principle behind the debate over "elasticity-functioning" versus "substrate porosity making sense" ^{[131,132,133,134]}. Beyond network porosity, our research has shown that the fibrous morphology between crosslinks—straight or crimped—also substantially affects the collagen macroscopic elasticity ^[87]. By distinguishing between viscoelasticity and poroelasticity as reviewed in Section 4, the viscoelasticity effect on cell behavior may well be due to the porosity of substrates. To bridge this gap, we suggest a framework that directly exploits the experimentally quantified fibrous matrix using imaging and graphical technology for the vicinal substrate around and contacting the cell, as was carried out in ^[135], with the addition of a continuum media model to represent the boundary conditions. Modeling the mechanical properties of soft biomaterials to dig into the micro physical fundamentals would propel advances in mechanobiology.

---

Acknowledgement

This study is funded by Grant-in-Aid for Scientific Research (C) from the Japan Society for the Promotion of Science (JSPS) (17K01352).

---

Conflict of Interest

The authors declare that there is no conflict of interest regarding the publication of this manuscript.

---