A new analytical method for estimating lumped parameter constants of linear viscoelastic models from strain rate tests

We introduce a new function, the apparent elastic modulus strain-rate spectrum, Eapp(ε˙)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$E_{\mathrm{app}} ( \dot{\varepsilon} )$\end{document}, for the derivation of lumped parameter constants for Generalized Maxwell (GM) linear viscoelastic models from stress-strain data obtained at various compressive strain rates (ε˙\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\dot{\varepsilon}$\end{document}). The Eapp(ε˙)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$E_{\mathrm{app}} ( \dot{\varepsilon} )$\end{document} function was derived using the tangent modulus function obtained from the GM model stress-strain response to a constant ε˙\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\dot{\varepsilon}$\end{document} input. Material viscoelastic parameters can be rapidly derived by fitting experimental Eapp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$E_{\mathrm{app}}$\end{document} data obtained at different strain rates to the Eapp(ε˙)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$E_{\mathrm{app}} ( \dot{\varepsilon} )$\end{document} function. This single-curve fitting returns similar viscoelastic constants as the original epsilon dot method based on a multi-curve global fitting procedure with shared parameters. Its low computational cost permits quick and robust identification of viscoelastic constants even when a large number of strain rates or replicates per strain rate are considered. This method is particularly suited for the analysis of bulk compression and nano-indentation data of soft (bio)materials.


Introduction
To date, different testing and analysis methods are used to derive quantitative mechanical properties for describing intrinsic material behavior and predicting material responses under specific loading conditions (Mattei and Ahluwalia 2016). Although mechanical characterization is often limited to the derivation of elastic properties, many materials, such as soft tissues and hydrogels, are known to exhibit viscoelastic behavior with time-dependent responses (Fung 1993;Mattei et al. 2014Mattei et al. , 2015bMattei et al. , 2017aOyen 2013;Tirella et al. 2014;Zhao et al. 2010). Creep and stress-relaxation tests are among the most common methods used to characterize material viscoelastic properties (Chaudhuri et al. 2015;Chin et al. 2011;Higgs and Ross-Murphy 1990). If a material is subjected to strain or stress inputs small enough so that its rheological functions do not depend on the input level, the material response is said to be within the linear viscoelastic range (LVR).
Another popular method to characterize material viscoelastic behavior is through dynamic mechanical analysis (DMA), based on applying a small amplitude cyclic strain (or stress) input on a sample and measuring the resultant cyclic stress (or strain) response (Gabler et al. 2009;Kiss et al. 2004;Soong et al. 2006). For a given sinusoidal strain input the resulting stress response will be sinusoidal if the applied strain is small enough so that the tissue can be approximated as linearly viscoelastic. Viscoelastic material response is characterized by a phase lag (δ) between the strain input and the stress response, which is comprised between 0°(purely elastic) and 90°(purely viscous). The dynamic mechanical properties are quantified using the complex modulus (E * ) which can be thought of as an overall resistance to deformation under dynamic loading. The complex modulus is composed of the storage (E , elastic component) and the loss (E , viscous component) moduli, which are additive under the linear theory of viscoelasticity (E * = E +iE ) (Menard 2008).
However, it is well known that these viscoelastic testing methods may pose some limitations in the case of very compliant materials, such as soft tissues and hydrogels . Indeed, both step-response and dynamic mechanical tests require the establishment of an initial contact between the sample and the testing apparatus to "trigger" measurements. Starting from this contact point, a given stimulus (e.g. step of deformation or force in the case of stress relaxation or creep, respectively, or cyclic stress or strain in the case of DMA) is applied to the sample under testing, then the resultant response is measured. As a consequence, force-or strain-triggered methods are likely to cause significant pre-stress on very soft samples, altering their status. Moreover, they generally require quite long measurement trials, which may result in sample deterioration during tests . To address these issues, we have recently proposed the epsilon dot method (εM) for deriving material viscoelastic constants through rapid measurements performed at different strain rates within the LVR (Mattei et al. 2015b;Tirella et al. 2014). However, as it is based on a global multi-curve fitting with shared parameters of stress-time data acquired at different constant strain rates, the method becomes computationally expensive when a large number of strain rates and/or replicates per strain rate are considered. Here, we introduce the apparent elastic modulus strain-rate spectrum, E app (ε), as a function to derive lumped parameter constants for Generalized Maxwell (GM) linear viscoelastic models from a single-curve fitting of experimental E app data obtained at differentε.

Methods
The most general form of linear viscoelastic lumped model (i.e. with constant viscoelastic parameters) is the GM model ( Fig. 1) and consists of a pure spring (E 0 ) in parallel with N Fig. 1 Representation of the GM model Maxwell arms (i.e. a spring E i in series with a dashpot η i ) thus defining a set of N different characteristic relaxation times (i.e. τ i = η i /E i ) (Roylance 2001).
The stress-time response of a generalized Maxwell model to a constant strain-rate input with amplitudeε is given by where N is the number of spring-dashpot Maxwell elements in parallel to the spring E 0 .
Since ∂ε/∂t =ε is constant, t = ε/ε can be substituted in Eq. (1) obtaining the stressstrain response: The tangent modulus (E t ), defined as the slope of the stress-strain curve at any specified stress or strain point (Zhang et al. 2002), is given by which is a function of both strain (ε) and strain rate (ε) or, dually, a function of time t = ε/ε. For very high (ε → ∞) and low (ε → 0) strain rates, the E t becomes independent of the strain and is, respectively, equal to the sum of all the spring constants present in the lumped model (namely the instantaneous modulus, or to E 0 (generally noted as the equilibrium modulus E eq = E 0 ) forε → 0 as expected from classical viscoelastic theory (Lakes 2009).
Hence, from Eqs. (2) and (3), linear viscoelastic materials-that is materials which can be represented by linear viscoelastic models combining lumped parameters (i.e. springs and dashpots) whose elastic and viscous constants do not depend e.g. on strain or strain rateexhibit non-linear concave-down stress-strain (or dually, stress-time) responses to a constant strain-rate input. These non-linear responses are dependent on the strain rate (i.e. the stress at a given strain increases withε due to viscous effects, Fig. 2) and their tangent (E t ) decreases non-linearly over strain (or time) from E inst (at ε, t = 0) to E eq (for ε, t → ∞) because of the time-dependency owed to viscous relaxation phenomena. Thus, the changes in E t are dependent on both the strain rate and, of course, the viscoelastic characteristics of the material tested. Zhang et al. showed that: (i) the tangent modulus for constant strain-rate tests, E t (t), is equivalent to the relaxation modulus from stress-relaxation tests, E(t) = σ (t)/ε 0 , where σ (t) is the stress measured over time in response to a step strain input ε 0 , and (ii) the plot of σ (t)/ε versus ε(t)/ε, whose tangent returns E t (t) = E(ε(t)/ε), is no longer dependent on the strain rate (i.e. curves obtained at different strain rates can be merged into one) for linear viscoelasticity (Zhang et al. 2002). Note that if σ (t)/ε versus ε(t)/ε curves collected at different strain rates cannot be merged into one, the viscoelasticity of the material under testing is non-linear and the linear superposition principle cannot be applied.
If a material is known to exhibit pure linear viscoelastic behavior, then in principle the lumped constants of a given linear viscoelastic model (i.e. spring E i and dashpot η i values) can be derived by fitting the non-linear stress-strain data acquired at a single constant strain-rate to Eq. (2). The strain-range considered in the fitting should extend up to the region in which the stress-strain tangent (E t ) becomes independent of the strain in order to ensure that all the viscoelastic relaxation phenomena have occurred and hence are included in the fitted dataset. However, many real materials are non-linearly elastic, viscoelastic or inelastic. Thus, they generally exhibit either strain-hardening (concave up) or strain-softening (concave-down) stress-strain responses to a constant strain rate. Most biomedical materials, including soft tissues and hydrogels, exhibit concave up σ − ε curves (Fung 1993;Jaspers et al. 2014;De Pascalis et al. 2014;Rimmer 2011;Vena et al. 2006), a clear indication of non-linear strain-dependent behavior. As a consequence, linear viscoelastic models cannot be used to describe the material mechanics. Conversely, in the case of concave-down stress-strain curves it is almost impossible to determine a priori whether the non-linear behavior observed comes from linear viscoelastic phenomena or non-linear material characteristics.
As discussed above, viscoelastic materials do not exhibit an initial linear portion on their σ − ε curve to calculate a Young's modulus. Nonetheless, from an engineering point of view, they are often treated as linearly viscoelastic in the region of small deformations Tirella et al. 2014;Yan and Pochan 2010). Hence, it is common practice to identify an initial linear region (called the linear viscoelastic region, LVR) for deriving an apparent (i.e. strain-rate dependent) elastic modulus (E app ), which is independent of the magnitude of the applied strain or stress. The LVR can be considered as the first nearly linear tract of the non-linear stress-strain response exhibited by linear viscoelastic materials: within this region the material's molecular arrangements are close to equilibrium and the mechanical response is essentially a reflection of internal dynamic processes at the molecular level (Choi et al. 2015). Since material rheological functions do not depend on the input level within the LVR, linear viscoelastic models can be reasonably employed to ana-lyze experimental mechanical data. Depending on the testing method chosen to characterize material linear viscoelastic properties, mechanical data within LVR are typically acquired at: (i) multiple stress or strain step-input amplitudes (in the case of creep-compliance or stress-relaxation tests, respectively), (ii) constant stress or strain rate (constant-rate tests like the epsilon-dot method,εM), or (iii) multiple frequencies (dynamic mechanical analysis, DMA).
There is no explicit mathematical form to define the E app since, by definition, a linear regression procedure is required to derive the E app from experimental σ − ε data obtained at constantε within the LVR. However, a good approximation of the E app can generally be obtained by evaluating the tangent modulus (defined in Eq. (3)) at ε * = LVR/2 strain. For instance, if the LVR extends up to ε = 0.05 then E app ∼ = ∂σ ∂ε | ε=0.025 (see Supplementary  Information SI 1). Under this assumption, the E app epsilon dot spectrum, E app (ε), can be derived from Eq. (3), obtaining The E app (ε) can be used to identify lumped viscoelastic parameters (i.e. E i , η i ) from experimental E app data obtained at different constantε as shown in Sect. 3.

Results
To investigate whether the E app (ε) fitting method proposed in this study returns viscoelastic constants similar to the originalεM described in Tirella et al. (2014), it is essential to maintain all sample, testing and analysis variables unaltered: only then can results be meaningfully compared (Mattei and Ahluwalia 2016). Thus, the same apparent elastic moduli reported in Tirella et al. (2014) for both 5% w/v Type A gelatin and PDMS (10:1 base to catalyst) tested in unconfined compression at different constant strain rates (n = 3 replicates per strain rate) were used in this work to obtain two experimental E app vs.ε datasets to analyze with the new E app (ε) fitting method. Moreover, to exclude any analysis-related variability factors other than the analysis method, the same LVR range and linear viscoelastic model previously used in Tirella et al. (2014) were used also for the E app (ε) fitting. Therefore, the E app vs.ε datasets obtained from Tirella et al. (2014) were fitted to the E app (ε) function of a Maxwell Standard Linear Solid (SLS) model (Eq. (5), obtained by substituting N = 1 in Eq. (4)), fixing ε * = LVR/2 = 0.005.
To avoid local minima in the fitting procedure, different sets of initial guess parameter values were used for the SLS viscoelastic constants to estimate (i.e. E 0 , E 1 and η 1 ). In particular, two different starting sets were defined respectively for gelatin (E 0 = 1 kPa; E 1 = 10 kPa; η 1 = 10 kPa s) and PDMS (E 0 = 100 kPa; E 1 = 10 kPa; η 1 = 100 kPa s) samples. Then an annealing scheme-multiplying and dividing each initial guess parameter value by 10-was adopted, leading to a total of 7 different initial guess sets. The minimum parameter values were constrained to zero to prevent the fitting procedure returning negative estimated viscoelastic coefficients. Comparisons between viscoelastic parameter values were made using the Student t-test, considering n = 3 (i.e. the number of replicates per strain rate) and setting significance at p < 0.05. The fitting procedure and the statistical analysis were implemented   Figure 3 shows experimental E app data and fitted curves for gelatin and PDMS samples tested in unconfined compression at differentε. The fitting results for Maxwell SLS models are summarized in Table 1, where E inst and E eq represent the instantaneous (i.e. E 0 + E 1 ) and equilibrium (E 0 ) moduli, respectively, while τ is the characteristic relaxation time, equal to η 1 /E 1 . To compare the computational cost between the new E app (ε) fitting approach proposed in this study (based on a singlecurve fit) and the originalεM (based on multi-curve global fitting with shared parameters), the computation time is also shown in Table 1.
Fitting results were independent of the set of initial guess parameter values both for the E app (ε) and theεM method ( Supplementary Information SI 2). No statically significant differences were found between the viscoelastic parameters estimated with the E app (ε) fitting and theεM for either gelatin or PDMS samples, demonstrating the equivalence of the two analysis methods (Table 1). Moreover, the significantly lower computational cost highlights the efficiency of the E app (ε) fitting method with respect to the originalεM.

Discussion
To overcome the issues of pre-stress and long experimental testing trials for viscoelastic characterization of soft materials using classical testing methods, we developed the epsilon dot method (εM) which is based on a series of rapid strain rate measurements (Mattei et al. 2015b;Tirella et al. 2014). However, theεM derivation of material viscoelastic constants involves a global multi-curve fitting of stress-time data obtained at different constant strain rates, performing chi-square minimization in a combined parameter space where viscoelastic constants to estimate are shared between the globally fitted curves . As a result, this method is computationally expensive and difficult to automate when a large number of strain rates and/or replicates per strain rate are considered.
Starting from the strain-rate dependent apparent elastic modulus (E app ) typically used for describing material viscoelastic behavior, we introduce the concept of the E app epsilon dot spectrum, E app (ε). This function enables the calculation of viscoelastic constants for Generalized Maxwell (GM) models through a single-curve fit of E app data obtained at different constant strain rate (ε) versusε, according to Eq. (4). The E app (ε) can be interpreted as the strain-rate domain equivalent of the storage modulus frequency spectrum, E (f ), typical of DMA analysis.
The E app (ε) method proposed in this study was validated versus the originalεM described in Tirella et al. (2014). As discussed in Sect. 2, unlike purely elastic materials, viscoelastic materials do not exhibit an initial linear portion on their stress-strain curve. Thus, the selection of the LVR is critical for any study dealing with viscoelastic material characterization and it is likely to affect the derived mechanical properties. According to our previous work on theεM , we here identified the LVR as the region in which stress varies linearly with strain with R 2 > 0.999. To obtain meaningfully comparable results, two experimental E app vs.ε datasets to analyze with the new E app (ε) fitting method were obtained from the same apparent elastic moduli reported in Tirella et al. (2014) for gelatin and PDMS samples. For the same reason, the E app (ε) fitting was performed considering the same LVR range (0 ÷ 0.01 strain range) and linear viscoelastic model (Maxwell SLS) previously used to derive lumped viscoelastic constants with theεM.
Viscoelastic parameters for a Maxwell SLS model estimated using the E app (ε) equation were found to be statistically similar to those obtained with theεM, demonstrating the equivalence of the two fitting approaches. Notably, the standard error of parameter estimations obtained with the E app (ε) fitting are slightly higher than those of theεM, because fewer data are considered in the fitting procedure. In fact, for each experimental measurement made at a given strain rate, the former method considers only n = 1 representative point (i.e. E app ) instead of the entire stress-time curve within LVR considered by theεM.
It should be noted that, although based on a non-linear function (Eq. (4)), the E app (ε) method can only be used to characterize linear viscoelastic properties. Indeed, the E app (ε) function is non-linear only with respect to the testing input variable (i.e. the strain rate,ε), but not to its viscoelastic parameters (i.e. spring E i and dashpot η i values), which are considered constants. This is similar to the function used to fit stress-time data in the original εM or to those used to fit storage and loss moduli obtained from DMA, which are both non-linear in their respective testing input variable (i.e. strain rate in the case ofεM or frequency in the case of DMA), but constant in their viscoelastic parameters. Therefore, likė εM and DMA, this method can only be used to analyze LVR data and describe material linear viscoelastic properties.
In conclusion, the E app (ε) fitting method presented in this work is robust and accurate and, more importantly, because a single curve is used to describe material behavior over the strain-rate spectrum, it is computationally more efficient than the originalεM (as demonstrated by the significantly reduced computation time, Table 1). Hence it could facilitate the adoption of strain-rate mechanical testing methods and subsequent calculation of viscoelastic parameters among researchers interested in characterizing viscoelastic (bio)materials. Nano-indentation methods are particularly suited for the automated derivation of material viscoelastic properties with this approach, since they allow the acquisition of E app data at several strain rates in different locations of the same sample (Mattei et al. 2015b). Besides nano-indentation, it is worth noting that the E app (ε) fitting method presented in this work can be adopted to analyze LVR stress-strain data obtained from any type of constant strainrate testing approach, regardless of the measurement length-scale, such as bulk compression (as in this work) or tension and macro-indentation. Potential applications of this viscoelastic characterization method range from mechanical and structural design for (bio)material and civil engineering (Jelen et al. 2013;Öchsner and Altenbach 2015), to tissue engineering (Mattei et al. 2017b;Tirella et al. 2015) and cell mechano-biology (Mammoto et al. 2013;Mattei et al. 2015a).