Comparative Analysis of Nonlinear Viscoelastic Models Across Common Biomechanical Experiments

Biomechanical modeling has a wide range of applications in the medical field, including in diagnosis, treatment planning and tissue engineering. The key to these predictive models are appropriate constitutive equations that can capture the stress-strain response of materials. While most applications rely on hyperelastic formulations, experimental evidence of viscoelastic responses in tissues and new numerical techniques has spurred the development of new viscoelastic models. Classical as well as fractional viscoelastic formulations have been proposed, but it is often difficult from the practitioner perspective to identify appropriate model forms. In this study, a systematic examination of classical and fractional nonlinear isotropic viscoelastic models is presented (consider six primary forms). Consideration is given for common testing paradigms, including varying strain or stress loading and dynamic conditions. Models are evaluated across model parameter spaces to assess the range of behaviors exhibited in these different forms across all tests. Similarity metrics are introduced to compare thousands of models, with exemplars for each type of model presented to illustrate the response and behavior of different model variants. The parameter analysis does not only identify how the models can be tailored, but also informs on the model complexity and fidelity. These results illustrate where these common models yield physical and non-physical behavior across a wide range of tests, and provide key insights for deciding on the appropriate viscoelastic modeling formulations.

Illustration of the storage and loss modulus for common biomaterials. A) Storage modulus for common biomaterials adapted from Budday et al. [21]. B) Loss modulus for brain [13,22,39], lung [118], breast [133], liver [88,122], skeletal muscle [141], myocardium [106,121], skin [62,79], artery [81,129], cartilage [42] and bone [18]. C) Examples of tissues in rheological testing -brain [21] in shear testing device, liver [9] in torsion device, myocardium in triaxial shear device [136], artery [24], skin in multiaxial device [52] and bone [70] 1 Introduction Accurate biomechanical model of tissues has a wide range of applications in medicine, such as patient-specific modeling, treatment planning, tissue engineering and in vivo tissue assessment. In particular, patient-specific models have been employed to understand a range of conditions, including mechanical properties and functions of the heart [6,7,25,59,83,112], mechanical function of cardiac valves [2,126], response to diseased abdominal aortic aneurysms and dissection [50,119,132], and collision mechanics in the brain [53,156]. Appropriate tissue models enable pre-operative planning through simulations and have been deployed for bone implants [29,55,107], artery stent insertion [5,40,54], and surgical training [93], amongst other applications. They are also critical for diagnostic methods such as magnectic resonance elastography (MRE) [82,94,133,134,148], which rely on models to translate wave motion into local tissue properties. In addition, characterization of tissues using such models also plays an important role in tissue engineering and fabrication of biocompatible materials [77,158]. A common thread through all applications is the need for accurate biomechanical models that encapsulate the stress-strain response over a range of loads and dynamic time-scales.
Development of appropriate biomechanical models involves replicating experimental results in a manner that elucidates the link between the stresses and strains. Tissues exhibit a broad range of mechanical responses (see Fig. 1), and have resulted in the development of various experimental approaches to mechanical testing. Tensile tests are some of the earliest and most straightforward experiments to be performed, with deformations being imposed on tissue samples either uniaxially (most commonly) or biaxially [33,60,87]. Compression tests have also been employed, with compression being applied either uniformly on a sample [103], or through local indentation [80]. Simple shear [37] and torsional shear [9,140] provide yet another mode of probing material response. Pure shear is the homogeneous non-rotational flattening of a body, in which the principal strains remain parallel to their respective principal stress axes during deformation. For an incompressible material, this can be achieved by holding one side at constant length while stretching or compressing one of the two remaining sides. For small deformations, simple shear and pure shear differ by rotation [72], but this does not hold true for large deformations [104]. Beyond the mode of deformation, experiments can utilize different dynamic loading/unloading protocols, illustrating a tissue's response to frequency [74,115] or relaxation behavior [44,150]. Alternatively, a sample can be subjected to deformation under a constant load, which characterizes the creep behavior of a tissue [80,116]. The observed behavior under these tests may be post-processed to describe the tissues response, informing the selection of an appropriate constitutive model.
Hyperelasticity has been widely investigated in biomechanical modeling for a range of tissues including but not limited to artery [45,51,66,68], ligament and tendon [69,131,153], valve [85,96,113,126], brain [20,75,98,102], myocardium [8,48,64,151], liver [28,47,123], and skin [3,14,57]. Since inception, the field has burgeoned with many different constitutive formulations including polynomial [95,99], exponential [34,43,147], logarithmic [68,139] and Ogden [109] forms (or some combination [28,145]). A review of different hyperelastic constitutive model forms in biomechanics can be found in [67,154]. However, despite evidence of viscoelasticity being observed since the early biomechanical works [43,44,87], it has been neglected in favor of hyperelastic only approaches due to focuses on quasi-static analysis, limitations in available experimental data, and computational considerations. While hyperelastic formulations have enabled studies of biomaterials, viscoelastic effects such as stress relaxation, creep, hysteresis, and variable frequency response are often ignored (Fig. 2). Viscoelastic modeling was first investigated by bringing together linear rheological elements that exhibit purely elastic or purely viscous behavior -springs and dashpots [87]. The simplest and earliest forms involve these two elements arranged in series -Maxwell model, or in parallel -Kelvin-Voigt model [27,30,38,56,100]. To match physiological responses, these models are extended with an increased number of rheological elements, e.g., the standard linear (Zener) model or the generalized Maxwell model [124,155,159]. Extension of these linear models to nonlinear viscoelasticity has also been presented in the literature, following varying generalizations [11,31,58,63,65]. These viscoelastic formulations provide a foundation for accurately capturing the rate-dependent behaviors of tissues.
An alternative to classical approaches are the class of fractional viscoelastic models (see [91] (1D) and [41] (3D)). Here, the first order derivatives common to classical models are replaced with fractional time derivatives, enabling the rate-dependency of the model to transition between elastic and viscous. This change in differential operator also introduces a relaxation spectrum with rate responses reflective of the power-law behavior observed in tissues [61,84,90,137]. Although this can be replicated using multiple branches in the generalized Maxwell model [155,160], the fractional approach achieves this using a single fractional derivative order α. While these models pose more complexity in their use for numerical simulations, efficient numerical approaches have recently been introduced that enable use of fractional models for tissue modeling [160]. The potential for these models to effectively simulate rate response in tissues and the development effective computational tools make fractional models another practical resource for describing tissue behavior.
In this work, we study a range of incompressible, isotropic and viscoelastic models across various testing protocols. The response of classical models (Maxwell, Kelvin-Voigt, Zener, generalized Maxwell) and their viscoelastic adaptions is investigated under oscillatory and long-term response (relaxation and creep) in uniaxial and biaxial stretch, compression, pure shear, as well as simple and torsional shear. We perform a methodical analysis on the impact of parameters on model behavior, with the aim of gaining an understanding on how these models can be employed and used to reproduce experimental data across a broad range of test conditions. This investigation enables a close look at these various viscoelastic forms and comparative analysis of their behaviors, providing a clear picture of how these models can be used to improve the modeling of tissue response.
This work is structured into several sections. Firstly, a brief review of fractional derivatives and kinematics is presented (Sect.

Methods
In order to describe the behavior of various models across testing protocols, several technical aspects need to be clarified. Firstly, fractional derivatives and kinematics are briefly reviewed in Sect. 2.1-2.2. The viscoelastic modeling approach followed in this paper is described throughout Sect. 2.3, while the testing protocols are kinematically defined in Sect. 2.4. The parameter variations considered for each model across testing protocols is systematically presented in Sect. 2.4.3. Lastly, details on the numerical implementation of fractional derivatives and stress computation are explained in Sect. 2.5.

Review of Fractional Derivatives
We start by reviewing the concept of fractional derivatives, which is important to establishing the fractional viscoelastic models in later sections. Fractional calculus started with the generalization of integral and differential operators from the set of integers to the set of real numbers (for reviews, see, Ross [125] and Machado [89]). Many different definitions for fractional differential and integral operators of arbitrary order have been introduced [117], with the Caputo definition being a popular choice for physics-based applications. For a ntimes differentiable function f and 0 < α ≤ n, the Caputo derivative can be written as, where α denotes the ceiling of α and is the Gamma function. For this work, we will only consider 0 ≤ α ≤ 1. The Caputo form (which can be related to other formulations like the Riemann-Louiville fractional derivative) has the advantage that The fractional derivative (for α / ∈ N), unlike integer derivatives, is a convolution of the past behavior of the function, making it particularly useful for incorporating the history of the stress of the material.

Review of Kinematics
To briefly review kinematics and kinetics in nonlinear mechanics [16,63,108,130], let Ω 0 ⊂ R 3 denote the reference configuration of a three-dimensional solid and Ω(t) ⊂ R 3 denote its configuration at time t . At every time t , it is assumed that there exists a material displacement, u : Ω 0 × [0, T ] → R 3 , relating the spatial position x ∈ Ω(t) to the material position X ∈ Ω 0 by where X i and x i denote the components of X and x. We also assume that, due to the incompressible nature of the tissue, there exists a hydrostatic pressure p : Ω 0 × [0, T ] → R 3 resisting volumetric change. This mapping in Eq. (3) is assumed to be bijective. The local deformation gradient F defining the mapping of Ω 0 to Ω(t) is given by where I is the identity tensor. The determinant J = det F defines the relative change in volume. The material stretch is characterized using C and B, denoting the right and left Cauchy-Green tensors [63,78], respectively, e.g., For isotropic hyperelastic biomechanical materials, the strain-energy function, = (I C , I I C , I I I C ), is often defined using three invariants of C [15], 1 I C = C : I, II C = C : C, I I I C = det C.
From the second law of thermodynamics, the strain-energy function can be related to the second Piola-Kirchhoff stress tensor (PK2) as which, in turn, is related to the first Piola-Kirchhoff, P, and the Cauchy stress tensors σ by push forward operations on S, In what follows, we focus on constitutive forms of the second Piola-Kirchhoff stress, S.

Constitutive Models of Viscoelastic Materials
In constitutive modeling of soft biological tissues it is common to express the hyperelastic stress in a material as the sum of the stress of its individual components (i.e., S e = S 1 + S 2 ...). These may be stresses due to anisotropic structures [64,66,98], hydrostatic effects [51,101,138], or varying degrees of nonlinearity [28,95,99]. Similarly, it is also common to express the viscoelastic stress, S ve , as the sum of the elastic and viscous components. Letting S e and S Q denote some symmetric tensor form, then a simple viscoelastic model could be S e and S Q are to be determined from the constitutive modeling of the tissue mechanical behavior and may or may not have the same form. While this encompasses models where elastic and viscous components are considered additive, a more generalized form can be written as Here a's and b's denote specific model parameters that enable consideration of a wide range of models (discussed further below, see Fig. 4). In this case, the viscoelastic stress is defined by a differential equation that must be solved (or approximated) to determine the current value of S ve , and is related to the elastic and viscous terms. Note that, the hydrostatic stress here is additive, and not encapsulated within the differential equation for S ve . This assumes that hydrostatic effects respond instantaneously without memory. This form encapsulates many of the classical viscoelastic formulations. However, extending to greater generality, the differential operators can be replaced with fractional variants, giving the general form, where D αa t denotes the Caputo derivative introduced in Eq. (1). Our main focus is to investigate different forms of S ve in both classical (α a , α b = 1) and fractional approaches (0 < α a , α b < 1) and examine their viscoelastic behavior over a range of kinematic conditions. There are many variables to be considered: the constitutive model for the elastic component, S e ; the constitutive model form for the viscous component, S Q ; the non-linearity of the hyperelastic and viscoelastic component; and the weighting of terms based on different values of a's and b's. 2 The variations considered, along with theoretical comparisons, are outlined below and summarized in Table 1.

Constitutive Form of S e and S Q
For the comparative analysis across viscoelastic models, we select a simple isotropic, exponential Fung-type model. Fung-type models are widely used in tissue modeling and their exponential form are versatile in capturing the behavior of various tissues in mechanical testing [44]. Here, we consider the specific form, S fung (c 0 , c 1 , C), characterized by two parameters c 0 , c 1 ∈ R + : S fung (c 0 , c 1 , C) = c 0 exp (c 1 (I I C − 3)) C. (9) In this stress definition, the parameter c 0 scales the stress-strain response, while c 1 enables varying degrees of nonlinearity. S fung can be directly derived from the strain energy function = c 0 2c 1 exp (c 1 (I I C − 3)) and was chosen based on previous experience. Other formulations could be considered, e.g., Mooney-Rivlin or Ogden-based forms as in Capilnasiu et al. [23]. However, the form used in Eq. (9) provides a good balance between model complexity and versatility. Further, we note that S fung maintains symmetry and objectivity, which is particularly important when considering differential forms as in Eq. (8). Using this constitutive form, we set with dev denoting the deviatoric operator, e.g. dev[A] = A − 1 3 (A : C)C −1 . 3 S Q was chosen to have the same form as S e to make it easier to interpret the extensive amount of results presented in the following sections. Given the generalized form shown in Eq. (8), we reduce 2 Note for practical purposes, parameters a 0 and b 0 are set to 1, leading to a 1 and b 1 setting the relative contribution of D α a t S ve and D α b t S Q , respectively. 3 Note that the deviatoric form of S Q was used for convenience to ensure that S Q = 0 when the body is unloaded and when prior strain history cannot be established at the start of the simulations.

Zener SLS integral form
Lin redundancy by setting c e 0 = c v 0 = 1. Moreover, while in general c e i and c v i may be selected arbitrarily, to reduce space of potential models, we chose these parameters to be equivalent for elastic and viscous tensor components.

Classical Viscoelastic Model Examples
This section focuses on classical viscoelastic models that are often encountered in literature. These models are summarised in Table 1 and are separated into their formal linear and nonlinear forms. The nonlinear analogues provide a form corresponding to the same modeling approach as the linear models but do not necessarily equate in the linear viscoelastic limit. To maintain object reference, the Cauchy stress σ is replaced by the PK2 tensor, S ve , in the nonlinear analogue. In the classical form, elastic contributions are replaced by S e and viscous contributions by S Q . For the classical forms, α a , α b = 1 in Eq. (8) as shown in Eq. (7). The first classical form in this study is based on the Maxwell model, composed of a spring and dashpot in series (see Fig. 3 A). The nonlinear formulation of the model, given below, represents a simplified differential form for S ve where Recall that S Q is defined as in Eq. (10). The model can be tuned through three parameters: a 1 , b 1 , and c v 1 , where a 1 and b 1 scale linearly the time derivative of S ve and S Q , respectively, and c v 1 scales the exponential power in Eq. (9). The Maxwell model is particularly useful in modeling relaxation behavior, but is known to be unable to predict typical creep behavior [100]. When held under constant stress, materials usually exhibit creep that gradually slow to a stop. However, this model displays a linear strain increase (constant rate).
An alternate form is the arrangement in parallel of a spring and dashpot, known as the Kelvin-Voigt model (see Fig. 3 B). Here, S e and S Q are defined as in Eq. (10), with the exponential parameters c e,v 1 controlling the nonlinearity of the model. Parameter b 1 acts as linear scaling factor, leading to a total of three parameters to tune the model. The Kelvin-Voigt model has the advantage that it can model creep behavior, but inherently it cannot predict stress relaxation. Instantaneous drop in velocity will result in an instantaneous drop in stress, whereas an exponential decay is normally observed in tissues (see Fig. 2).
The standard linear solid model, also known as the Zener model, combines two springs and a dashpot. One possible arrangement is having a spring in series with a Kelvin-Voigt configuration or, as considered here, having a spring in parallel to a Maxwell configuration, as shown in Fig. 3 C. In this configuration, the nonlinear Zener model becomes As before, S ve , S e and ∂S Q /∂t are the viscoelastic, elastic and viscous PK2 stresses, respectively. Four parameters are used to describe the model. Firstly, c e,v 1 scale the exponential power and thus have a nonlinear influence on the model, while a 1 and b 1 act as linear scaling factors. By having a more complex arrangement of rheological elements, both in series and in parallel, the Zener model bypasses the issues encountered by the Maxwell and Kelvin-Voigt models in creep and relaxation, respectively.
Another development of classical rheological models is the generalized Maxwell model with n branches, as depicted in Fig. 3 D. Note that the Zener model is a particular form of the Generalized Maxwell model with only one branch. The benefit of having multiple branches is that they allow for a more diverse set of rate responses within the model, enabling a more tailored relaxation response [1,127,128]. For compactness, this model is presented below in its integral form, as A total of 2n + 2 parameters are needed in this case. Firstly, the exponential power scales c e,v 1 trigger a nonlinear response in the formulation of the PK2 stresses S e and S Q . Then, each branch is defined by the relaxation time τ k and relaxation weight b k . The generalized Maxwell model is the most comprehensive classical viscoelastic model, providing flexibility in representing decaying relaxation spectra. However, this generality comes at the added cost of an increased number of parameters.

Fractional Viscoelastic Models
Fractional models provide an alternative to classical models in modeling viscoelasticity. From a rheological standpoint, a fractional element can be imagined as a spring-pot, analogous to a more complex set of Maxwell elements arranged in parallel [160]. Thus, it is functionally similar to Eq. (14) with pre-established weights and decay constants that are not free parameters, but tied to the underlying fractional order. This is achieved using the fractional derivative (Eq. (1)) with the order 0 < α < 1. As in the case of the classical models, the fractional models discussed here are summarised in Table 1. There, the linear form of the models is provided, alongside the nonlinear analogues which are further presented in this section. The first fractional model considered here will be referred to as the fractional Kelvin-Voigt, which provides a direct alternative to the generalized Maxwell model. The nonlinear fractional Kelvin-Voigt model is described by Note the similitude between the fractional and classical Kelvin-Voigt models shown in Eqs. (15) and (12). If the fractional order is chosen such that α b = 1, then Eq. (12) is recovered, but the benefits of the fractional model are lost and relaxation issues arise. This is because, as α b increases, the relaxation spectrum of the model increasingly weights more recent events over the past. Overall, the fractional Kelvin-Voigt model captures the same complexity as the generalized Maxwell model (with preset assumptions about the distribution of the spectra) and uses four parameters -the scaling parameters b 1 , the exponential scaling powers c e,v 1 and the fractional order α b . Another fractional model analogue to a classical configuration is the fractional Maxwell model. This is derived by replacing the dash-pot with a spring-pot, leading to a series configuration as seen in Fig. 3 F. Following the form shown in the classical Maxwell (Eq. (11)), the fractional version replaces full derivatives with fractional counterparts, e.g., Five parameters define the model: the linear scaling factors a 1 and b 1 , the fractional orders α a and α b and the exponential power scaling c v 1 . For non-extreme α a , α b values (i.e., α a , α b < 1), the fractional Maxwell model may change in the creep description, but this response will be investigated later in this study.
Lastly, here we devise a comprehensive fractional model that combines the fractional Kelvin-Voigt and Maxwell formulations, with the liberty of having a different fractional order derivative of the total viscoelastic stress S ve and the viscoelastic stress S Q . In this case, which will be referred to as the fractional Zener model. Here, six parameters are needed for the model: the fractional orders α a and α b , the linear scaling terms a 1 , b 1 , and the exponential power parameters c e,v 1 . This is the most general model form, from which the other models can be obtained (see Fig. 4), e.g., the nonlinear fractional Maxwell by disregarding the S e contribution (i.e., setting b 0 to 0), or even the classical form by additionally setting α a = α b = 1. The transformation to the generalized Maxwell model is not immediately obvious. However, by recalling that the nonlinear fractional Kelvin-Voigt is an alternative approximation to the generalized Maxwell model, then setting a 1 = 0 provides the correspondence.

Common Biomechanical Testing Paradigms
To examine the capacity of the different classical and fractional viscoelastic models to exhibit realistic behaviors, a series of common biomechanical testing paradigms were considered. While there are numerous testing protocols in the literature, this work focused on strain loading cases, including uniaxial extension/compression, biaxial extension (with pure shear as a special case), simple shear, and torsion, as well as stress loading such as creep. To capture viscoelastic responses, different modes of cyclic (with varying frequency) or static loading were also considered.

Dynamics of Strain and Stress Loading
The impact of viscoelasticity on the mechanical function of biological tissues and organs can be divided into the short and long-term responses. In vivo periods of short-term cyclic loading (e.g., pulsation due to flow, breathing, walking, etc) and long-term static loading (e.g., diastolic blood pressure, standing, sitting, sleeping, etc) can be observed. As a consequence, experimental rheological studies have been designed to probe a broad range of tissue responses. Capturing both cyclic and static responses, we introduce three primary dynamic functions, covering common dynamics used in rheological experiments, a step (t, τ P ) = max 0, min 1, denoting sinusoidal temporal variations, sawtooth variations and periodic step functions.
Here, τ P denotes the period of oscillation, which can be extended for long-term static loading. Note that dynamic functions are bounded in time, e.g., |a k (t)| ∈ [0, 1], ∀t ∈ R. Further, a step is introduced whereby 10% of the period is spent ramping up or down the load, 40% of the period spent statically holding the load, and 40% of the period spent recording recovery.
In this paper, we consider both strain loading and stress loading. Focusing initially on strain loading, we consider both stretch and shear. Using the dynamic function of Eq. (18), the stretch, λ, and shear, γ , can be defined as, Here λ max denotes the maximal stretch or shear achieved. Stress loading is also applied in the literature, whereby a preset force is prescribed and the material stretch or shear is recorded. In this paper, we focus on creep response, in which a certain set load, T , is applied to the material and released, i.e., Here τ P denotes the period of on/off loading and unloading, and T max denotes the total load applied. In this work, only step loading is considered, reflective of the typical dynamics applied in the literature. The following sections outline the kinematics applied in strain loading. In this case, the stretch, λ, and shear, γ are referred to without parameters or superscripts for ease.

Kinematics of Simple Deformations
This work focuses on the viscoelastic response of models to the 4 pure deformation kinematics, uniaxial, biaxial, simple shear, and torsion. The isotropic nature of the mechanical response considered removes the need for consideration of different loading directions, reducing the number of kinematic variants considered. For uniaxial deformation, we will examine both the tensile and compressive cases, where the deformation gradient in vector notation is given by where λ 1 denotes the applied stretch or compression (applied in theê 1 direction). The biaxial deformation will be restricted to a general case, where λ 2 > λ 1 , and pure shear. Biaxial deformation in general is given by where one definition of pure shear can be expressed as which is more commonly found in experimental studies as it is much simpler to be applied mechanically [104]. For testing the mechanical response under shear, we consider mainly the simple shear case given by Torsion testing is a common method for probing the viscoelastic properties of 3-dimension tissues (where the membrane assumption is not applicable). Here a block or more commonly a disk of tissue is placed between a stationary plate and a linear actuator. An initial small compressive load is first applied and then the actuator will twist cyclically at a predetermined frequency. The deformation applied consists of both compressive and extensional loading as well as shearing, and the deformation applied is also heterogenous, as it shears more away from the central axis of the tissue. Simulating the torsion testing is thus a good example application to examine the mechanical response of the viscoelastic model forms.
We do note that the method of attachment is a complex issue and the deformation exhibited is often non-ideal (i.e., the deformed tissue will bulge at the middle). Our goal is not to replicate the full experimental setup and will as such assume ideal deformations, i.e., the geometry of the tissue will remain perfectly cylindrical. The deformation gradient for such a setup is given by where γ is the shear applied by the twisting angle , R is the initial radius of the specimen, and λ z is the initial compressive loading applied.

Classical and Fractional Viscoelastic Model Analysis
In order to understand the envelope of behavior for each model, simulations were performed across all testing paradigms using varying model parameters. As previously mentioned, parameters c e 0 and c v 0 were set to 1, to avoid redundancy with the b 0 and b 1 parameters. The exponential parameters, c e 1 = c v 1 = c 1 , were set to be equal for simplicity. Without loss of generality, parameter a 0 is set to 1 (Eq. (8)). For analysis purposes, the elastic part parameter b 0 was also set to 1 (for the Zener and Kelvin-Voigt forms) or zero (for the Maxwell model forms). While changing b 0 alters the material response, it can be shown that the behavior can be replicated with a b 0 = 1, re-weighting of a 1 and b 1 , and a scaling of the total stress response (except in the case where b 0 = 0). In our analysis, we focus on comparing normalized stresses/strains (i.e., divided by their respective maximum value, which after normalization becomes 1), simplifying the need to consider all three parameters. Following the same reason, b 1 = 1 for all Maxwell model forms.
As a result, a 1 and b 1 become the principal scaling parameters responsible for the relative contribution of the viscoelastic components. Below we list the parameter variations considered in this analysis: Examining model behaviors, we aimed to address model responses for varying testing paradigms. Each testing paradigm includes the variation of kinematics, load level, cycling period, and cycling wave function. The kinematics include uniaxial, biaxial, pure shear, simple shear, torsion, uniaxial relaxation, and uniaxial creep. The loading parameter for each testing paradigm varies as following For the uniaxial case, a frequency sweep of the dynamic tensile modulus at small strains (1%) was also conducted. Here, frequencies of 10 −3 , 10 −2 , 10 −1 , 1, 10 1 , 10 2 , 10 3 , 5 × 10 −3 , 5 × 10 −2 , 5 × 10 −1 , 5, 5 × 10 1 , and 5 × 10 2 Hz were calculated. Loading was applied using the sinusoidal wave with λ max = 1.01 and the nonlinearity of the model was chosen to be c 1 = 0.5. Only uniaxial deformation is considered due to the limited range of strain and the nonlinearity of the model. For determining the storage and loss moduli over 10 −3 to 10 3 Hz, the stress and strain responses were fit to a sine function on the 5th oscillatory cycle. The Cauchy stress over time, σ , and strain over time, , were fit to using the optimization library from Mathematica [157] (where ω = 2π/τ P ). Following Meyers and Chawla [100], the storage and loss moduli were computed as Due to the significant amount of data available, we present the data analysis by focusing on outliers and the best case parameters for each model form. We devised an approach by computing a similarity metric between the models, examining how this metric varies with different model forms and parameters with a focus on the worst values. Given that all existing data was simulated, the sum of the residual squared is used to measure the difference between two curves. This value was then normalized by the mean total variance between the two curves. Thus, the metric is given by where f i is the curve of the normalized stress or strain component (by maximum positive value), respectively. For each kinematic, load magnitude, loading frequency, and material nonlinearity, we first compared how similar each pair of model forms is, to get an overview of model responses under the tested conditions. While this examination provides insight into contrasting or similar behaviors in models, it does not capture whether models exhibit physical behavior, such as those seen in Fig. 2. To get an impression about this, a representative response was selected from the Fractional Zener models that exhibited realistic responses under the different testing modes. All other models were then examined to find which best matched the exhibited response, providing a way of determining whether physically reasonable responses were viable for a given model. Given what was observed, we focused on presenting the result in terms of the overall similarity patterns and best case results for the kinematics, loading level, and frequency response.

Numerical Approximation of Model Responses
The above examples were implemented and simulated using Python3 (Python Software Foundation, https://www.python.org/). The strain loading tests were straightforward to implement using the NumPy library for support of array operations [110,144]. The fractional derivative were computed using the Prony series method [160]. To summarize, this approach is equivalent to using N Maxwell elements arranged in parallel, each with their own stiffness β k and time constant τ k (k = 1 . . . N), chosen to optimally approximate the fractional derivative for a given order. Supposing we wish to take the α-order fractional derivative (α ∈ [0, 1]) of the temporal tensor function A : [0, T ] → R 3×3 at any time t ∈ [0, T ], then the Prony series approximation to the true Caputo derivative in Eq. (1) can be written aŝ where Q k : [0, T ] → R 3×3 is an introduced intermediate variable corresponding to the kth Maxwell element. Discretization ofD α t was accomplished applying standard backward Euler. The resulting discrete approximation at t n = n t , denoted asD α n , can be written aŝ Here the computational cost stays constant over time and the storage of historical values is limited to the size of A n−1 and Q n−1 k (for each k = 1, . . . N). Applied to Eq. (8), the fractional operator must be applied to both S ve and S Q . For test cases when the temporal dynamics of S Q were known, we note that the backward Euler derivative in Eq. (30) was computed using a more accurate approximation (using central difference with h = max(λ, 1)( 1/3 ) where is the machine epsilon for double precision numbers). After isolating the correct terms, the stresses at the nth step are calculated using The heterogeneity in strain for the torsion example (Eq. (25)) requires the computation of the stresses over the top surface of which the traction is applied. Taking into consideration the observations from Mazzia et al. [97], the stresses were computed at quadrature point on a disk determined using the Gauss-Legendre rule in r (n = 8) and mid-point rule in θ (n = 8). For these series of tests, t was chosen to be 0.01 times the period of one loading cycle.
For the uniaxial creep simulations, a forward simulation was not possible due to the combination of applied tension, incorporation of strain history and solving of the differential equation Eq. (8). Instead, at each time step, an inverse problem had to be solved whereby the applied traction (see Eq. (20)), T e 1 , on the surface with normal e 1 , had to balance with the stress at an unknown stretch, λ 1 , e.g., Without a defined λ 1 (t), (8)) was also computed implicitly by backwards Euler. Solving for λ 1 was numerically unstable for some of the classical forms with larger time steps, where convergence was not consistently achievable with Newton-Raphson iterations. Taking into account the time and data storage due to the number of simulations being ran, dt was chosen to be 0.005 times the period of one cycle and a more rigorous 2 layered approach was taken for solving Eq. (32), with Brent's method [19] initially and a differential evolution step if Brent's method fails. For Brent's method, the λ n−1 1 was taken as one bound of the initial condition. The opposing bound was found by repeatedly subtracting the previous step size times the sign of the residual from λ n−1 1 until a point where a residual of opposite sign was found. Brent's method with tol = 10 −10 was ran for 100 iterations. 4

Results
In our analysis, 4, 216 viscoelastic models were compared and tested across 144 testing conditions. Due to the amount of data, selected results are presented in the main body of the paper highlighting key outcomes and conclusions. Additional results are presented in the Supplementary Materials. Further, modeling results were organized hierarchically into an ordered list so that the extensive amount of data can be examined systematically. The classical and fractional forms are then further divided into three groups, Maxwell, Kelvin-Voigt and Zener, with the Zener forms further divided into 8 groups corresponding to a 1 = 0.01 to 0.8. The 8 Maxwell forms are listed from a 1 = 0.01 to 0.8, the 6 Kelvin-Voigt forms are listed from b 1 = 0.2 to 20, and the 6 forms from each Zener group are listed from b 1 = 0.2 to 20 throughout the rest of the analysis. This order of model form is maintained throughout the presentation of the results.
To organize the results, the similarity scores between the model forms were reordered in this fashion and compiled into mosaic plots. To find the comparison between two specific models, the position of each model on the legend can be found by 1) first picking out whether the modeling approach is classical or fractional, and if fractional what parameters α a and α b were used, then 2) the specific form of Maxwell, Kelvin-Voigt, or Zener, and finally 3) the parameters a 1 and b 1 of the model (Fig. 5A). Their positions on the legend points to a specific block within the mosaic plot which is shown by the red and blue arrows in Fig. 5B. The color of each block represents the similarity metric (Eq. (28)) computed from the gray lines representing the difference between the blue and red curves, which corresponds to the color of the models selected in Fig. 5A (Fig. 5C). In the following sections, we examine how the responses of these model forms change in different kinematics testing setups -with the maximum strain applied, the nonlinearity of the material behavior, and the frequency of loading. Figure 6 shows how the different model forms compare to each other under sinusoidal cyclic uniaxial stretch of 1/1.2 to 1.2 and a nonlinearity parameter of c 1 = 2. Each cell of the array plot shows the similarity score (Eq. (28)) between two specific model forms. The diagonal indicate the comparison between the same forms (with a score of 0), and are thus colored black. The plot is symmetric by definition and thus the top triangle is used for the response at 5 Hz and the bottom triangle is used to show the response at 0.05 Hz (0.5 Hz were also tested but not shown). The results show significantly different behavior at each frequency.

Model Response Under Varying Kinematic Conditions
Notably, the classical forms show significantly higher scores on average. At 0.05 Hz, the classical Maxwell forms in particular show significant differences with the other forms (Fig. 6). The classical Kelvin-Voigt and Zener forms compare well with the fractional forms, except when b 1 is small. This is due to the fact that small b 1 values reduce to a pure elastic  28)) and compiled into mosaic plots in Fig. 6 and for other kinematics shown in the Appendix. A) Illustrates how to find the location for the comparison between two models by 1) looking up whether the modeling approach is classical or fractional and if fractional, then what parameters α a and α b were used, then 2) the specific form of Maxwell, Kelvin-Voigt, or Zener, and finally 3) the parameters a 1 and b 1 of the model. B) Once their position on the legend are found, they point to a specific block within the mosaic plot. C) The color of each block represents the similarity metric (Eq. (28)) computed from the gray lines representing the difference between the blue and red curves, which corresponds to the color of the models selected in A) response on the right hand side. The fractional forms are generally very consistent, especially for its Zener forms. On the contrary, some of the trends seem to reverse at 5 Hz. The classical Maxwell form behaves much more comparably with the other forms, while the Kelvin-Voigt forms and Zener forms with small a 1 values are outliers (note Zener forms with low a 1 values reduce to the Kelvin-Voigt form). The fractional forms also seem to bias against high b 1 values more often. These types of behaviors illustrate how the b 1 terms become more dominate at high frequencies and how the a 1 term can heavily dampen the rate dependent change in response. To examine more closely, we compared each form against a representative specimen (Sect. 2.4.3), plotted their similarity score (Fig. 7 left) and the best case normalized stress strain curves for the classical Maxwell, Kelvin-Voigt and Zener forms (Fig. 7 right).   Fig. 2, the red vertical axis is for the normalized stress, the blue horizontal axis is for normalized strain and the green alternative is for normalized time.
The most important difference in trend is the swap in behavior between varying the a 1 parameter and the b 1 parameter. At 0.05 Hz, each Kelvin-Voigt and Zener groups are arranged in a V-shaped pattern indicating an optimal response for a given b 1 value within each group. The Maxwell forms in this case tend to occur in a cluster with the trend of decreasing scores with a 1 . On the contrary, at 5 Hz, the Maxwell forms are arranged in the V-shape while the Kelvin-Voigt and Zener forms are clustered. Interestingly, the Kelvin-Voigt and Zener groups form a V-shape showing the optimality for a specific a 1 value. We can observe this more clearly by examining the shape of the stress-strain curves (Fig. 7). The classical Maxwell model at 0.05 Hz leads to the most dissimilar curves across all tests, reflecting the high scores observed. The hysteresis is larger than for the other models. Additionally, there is a negative slope around the 0% stretch, as opposed to the other models which yield a positive slope. In contrast, the observed classical Kelvin-Voigt (blue) curves switch place with the Maxwell (red) at 5 Hz, yielding an abnormal response with high hysteresis.
We next repeated this analysis for the remaining kinematics, and the results for 0.05 Hz and c 1 = 2 are shown in Fig. 8. Broader frequency investigation of each kinematics, as presented in the uniaxial case in Fig. 6 and 7, is presented in the Supplementary Materials. In  The abnormal behavior of the classical Maxwell model can be observed here for all kinematics. In the relaxation test, although the decay part of the curve shows a power-law decrease, the plateauing part of the stress quickly reaches 0, which is atypical of tissues (Fig. 2). For the unloading path, due to the viscous stresses and absence of elastic ones, the Maxwell model produces the largest resistance when the strain is reduced from λ = 1.2 to 1, leading to a stress undershoot before decaying again to 0. The fractional adaption of the Maxwell model can be tweaked to produce better results except for the creep test, where the dissimilarity would remain high. Under several cycles of applying a constant load and then unloading, the Maxwell models would lead to a continuous increase in strain, with no recovery.
The classical Kelvin-Voigt and Zener models show similar results across the tests' range, particularly in the oscillatory tests. Note, however, that for the uniaxial test using a saw-type loading, a jump in the stress level is recorded for the classical Kelvin-Voigt model. This is because the abrupt change between loading and unloading, under the derivative, leads to a sudden positive to negative stress change, which is only partly corrected by the elastic term. The classical Zener model provides a better continuity in this sense. For the relaxation test, both models decay to the elastic stress level, but with the Kelvin-Voigt form decaying instantaneously at the end of the loading phase. As the strain decreases, the viscous part in both models leads to an undershoot, before the stress settles to 0. Note that this is smaller than for the Maxwell model due to the presence of the elastic term. For the creep and release test, both classical Kelvin-Voigt and Zener models see a gradual strain increase during the constant load phase, then a sudden, almost linear decrease during the recovery. Over multiple creep and release tests the recovery strain drifts higher. A fractional approach of both Kelvin-Voigt and Zener models improves the similarity with the representative curves across all tests. Although pictured only for the fractional Zener (chosen as representative), the relaxation curve decays gradually towards the elastic limit, as exhibited in tissues in literature (Fig. 2). The strain decrease leads to a stress undershoot, which slowly relaxes towards 0. The creep and release test shows a gradual strain increase, followed by a gradual, non-complete recovery, which drifts higher across multiple tests. The low dissimilarity score suggests that the same can be inferred about the fractional Kelvin-Voigt model. Note that the creep test's similarity scores show a clearer differentiation of parameters that are suitable for reproducing a representative behavior. Therefore, this test can help single out optimal parameters from a wider range of satisfactory combinations.

Optimal Model Form for Material Nonlinearity
Here, we examine the models' behavior for two exponential parameters (c 1 = 1 and c 1 = 4), for two stretch ratios (λ = 1/1.1 to 1.1 and λ = 1/1.4 to 1.4). In Fig. 9, the impact of the exponential parameter and deformation level are shown for the uniaxial test, with the other tests being presented in the Supplementary Materials. The organisation of Fig. 9 is similar to before, with the left column showing the similarity metric of the classical and fractional models (for some parameter combinations) compared to a representative case, and the right column showing the normalised best classical models curves compared to the same representative curve. For both stretch ratios, an increase in the exponential parameter c 1 leads to noticeable differences. For a stretch of λ = 1.1, the models' response deviates more from an elliptical path. However, for a larger stretch (λ = 1.4), not only does the general trend of the curves change, but the hysteresis also diminishes, relative to the peak stress level achieved. Among the models, it is once again evident that the classical Maxwell formulation for low frequencies leads to responses that are not representative of those observed in tissues, Interestingly, the comparison of model responses is not significantly affected by the material nonlinearity (c 1 parameter). However, a difference in the optimal parameter values (a 1 and b 1 ) are observed at larger loading strain. This seems to be due, in part, to the larger strains being more capable at changing the effective nonlinearity of the material response than simply the exponent c 1 yet a fractional approach would improve this aspect. The classical Kelvin-Voigt and Zener models appear to approach well the representative fractional Zener curve for a stretch of 1.1. However, at large stretches, the classical forms undershoot when switching from loading to unloading. A fractional approach would help in this sense, as indicated by the similarity metric. Overall, here we infer that by increasing the deformation level within a test can improve the identification of the exponent parameter.

Optimal Model Form for Frequency of Loading
To better investigate the effect of the frequency (rate) of loading, we examined the storage and loss moduli in the range 10 −3 -10 3 Hz. The simulations were reduced to the uniaxial case with a maximum strain of 1% and a nonlinearity of c 1 = 0.5. Figures 10 and 11 show the storage and loss moduli (E and E , respectively) of the classical and fractional models, where the layout indicates the varying parameters. Note that, due to the models considered being incompressible isotropic, the trends describe the shear storage and loss moduli (G  For the remaining forms, the scaling is shown on the right. The Maxwell form curves correspond to different parameters a 1 , rather than b 1 as for the other model forms, and are thus plotted with dashed line. The order is similar to Fig. 6. From left to right, Maxwell, Kelvin-Voigt and Zener model forms with different a 1 are considered (either classical or fractional). For each individual plot, multiple parameters are shown considering the other free model parameter. Note that Maxwell plots use dashed lines to indicate that a different parameter is being varied as well as the axis range on the left side, while all other plots use the axis range on the right are located on the right hand size of the 6 th column. With this layout in mind, it can be observed that the fractional Kelvin-Voigt response is repetitive between row groups 2-3 and 4-5, as the varying parameter α a does not appear in the formulation of these models and thus does not influence the response.
Typical frequency behavior varies across tissues and range of frequency employed, but some general trends can be identified. Often, both storage and loss moduli increase with frequency following a power law [35,105], which may transition from a weaker to a stronger power law when the upper range of the frequency is increased [17,26,62]. Alternatively, both moduli may appear as increasing with frequency and then plateauing [86], and the loss modulus may even follow a descending trend after an initial increase [13,88,142]. Across the board, we can see that some models or parameter combinations do not yield any of the expected trends. A fractional order α a larger than α b leads to a decrease of the storage modulus towards the end tail of the frequency range, while the loss modulus vanishes after reaching a peak value. Even for cases when α a ≤ α b , when the linear scaling parameter b 1 is smaller than a 1 in the Zener models, the storage modulus follows a decreasing trend or a decrease followed by an increase, while the loss modulus becomes negligible for most of the frequency spectrum. On the same note, an increasing a 1 parameter in the Maxwell models (where α a = α b and b 1 = 1) leads to an increase and then decrease of the loss modulus even below the initial levels reached at low frequencies. Therefore, for models that employ derivatives on the viscoelastic stress as well as on the total stress, it is advisable that the order of the viscoelastic stress derivative is higher than that of the total stress. Bagley and Torvik [10] found that, for the fractional Zener model to be well-defined across any frequency range, α a and α b must be equal. It is also preferable that the scaling parameter of the viscoelastic stress be higher than the one of the total stress, especially if the fractional orders are equal. If these conditions are satisfied, then varying the a 1 parameter can help modulate the frequency at which the storage modulus plateaus or changes from a weak to a strong power law behavior, and the loss modulus plateaus or reverses from an increasing to a decreasing trend. The maximum storage and loss moduli achieved could be altered through the a 1 parameter, but parameter b 1 is more impactful in this sense. The Kelvin-Voigt models capture a power law behavior of the loss and storage moduli, which may be tailored to transition from a weaker to a stronger power law through the manipulation of the fractional order α b and the linear scaling b 1 in the case of the storage modulus. However, it does not appear to capture a plateauing storage modulus behavior, or an increase followed by a decrease for the loss modulus.

Discussion
This study analyzes the behavior of the nonlinear classical Maxwell, Kelvin-Voigt and Zener models and their viscoelastic adaptions under typical kinematic conditions imposed in rheology testing. An extensive parameter examination was carried out, with models relying upon up to 6 parameters, each spanning several values. Section 3 presented the models' behavior under specific conditions or parameters, and here we aim to overview and gather the knowledge gained from those analyses. We will start from the fractional Zener model, as it provides the most generality, and then consider the simplifications that lead to the classical, Maxwell or Kelvin-Voigt forms (Fig. 4), and their impact on the tests studied.
The most comprehensive model considered in this study is fractional Zener (Eq. (17)) and entails 6 parameters. Its versatility allows it to capture typical tissue behavior in oscillatory, constant stress, constant load and frequency tests. The similarity metric for the fractional Zener model indicates an extensive choice of parameters that can be employed to replicate typical tissue behavior. However, this is also an indicator on the extensive rheological testing that needs to be performed in order to identify a unique set of optimal parameters. Indeed, under limited frequency regimes, the fractional Zener model can be seen to behave similarly to simpler models (depending on the test considered). The classical Zener model effectively reduces the fractional model complexity by setting both derivative orders α a and α b to 1 (Eq. (13)), and thus relying upon 4 parameters. This simplification leads to significant changes in the relaxation and creep tests. During a constant applied strain, the model decays rapidly to a constant (elastic) stress level. The short relaxation time is atypical of tissues, as seen in Fig. 2. The creep shows a characteristic tissue response during the constant load phase but not during the release, where the strain decreases linearly. Therefore, the classical Zener model might be ill suited to model tissue behavior under high impact (e.g., abdominal organs in a car crash), or the heart for growth and remodeling problems, where loading triggers the growth.
A further model simplification would be achieved by dismissing the elastic stress and setting b 0 = 0, leading to a classical Maxwell model. In all the curves shown at a small loading rate (0.05 Hz in Figs. 8 and 9), the model failed to capture a representative behavior, mostly due to excessive hysteresis. Therefore, it is unsuitable to be employed with slow rate behavior, e.g., the liver response under respiratory motion. Nonetheless, at faster deformation rates, the hysteresis diminished (Fig. 7). The fractional approach of the Maxwell model also appears to be able to produce representative tissue behavior, as it arises from the similarity metric, and is not restricted by the deformation rate. However, the pitfall of both classical and fractional Maxwell models remains their inability to predict creep behavior, as the strain keeps increasing indefinitely and there is no strain recovery during the offload periods. Hence this class of models is to be avoided for scenarios involving persistent mechanical stresses.
Another avenue of simplifying the Zener model is by setting the a 1 parameter to 0, leading to the Kelvin-Voigt type of models. This is particularly attractive because the ordinary differential equation becomes easier to implement and solve. The classical form exhibits the expected limitations in modeling relaxation behavior, which appears to decay instantly. Moreover, at higher frequency rates in oscillatory tests, the hysteresis amplifies to non-physiological levels (Fig. 7). This also reflects in the frequency behavior over a wider range, as the storage modulus remains constant but the loss modulus has an increasing trend (Figs. 10 and 11). Therefore, the classical Kelvin-Voigt model should not be employed with fast rate phenomena (e.g., elastography applications). The fractional model approach is better suited for relaxation and frequency investigation for both slow and fast rates, and appears to be an attractive alternative to the fractional Zener model, from numerical and parameter fitting process perspectives. Its limitation arises from the frequency spectrum investigation, which reveals a continuous increasing path for both the storage and loss moduli. The frequency range over which these moduli increase according to a power-law, before plateauing, varies across tissues -e.g., 0.01-10 Hz for liver [88] and cardiac thin filaments [142], 0.1-100 Hz in kidney [105]. Therefore, the use of the fractional Kelvin-Voigt model may have limitations in some elastography applications, or for modeling cardiac tissue dynamics.
Rheological testing of tissues has grown over the decades. It is now acknowledged that tissues generally exhibit nonlinearity and viscoelasticity, among other more specific traits. When planning rheological testing, it is advantageous to have an insight into how different protocols might accentuate certain features and help towards a better identification of model parameters. In this study we saw that, for an isotropic material, uniaxial and shear deformation lead towards a similar behavior being exhibited (Fig. 8), hence choosing a single test out of the two could help avoiding redundancy. Biaxial stretching enables a better differentiation between the quasi-linear behavior and the transition to nonlinearity, hence being preferable for characterizing a model's nonlinearity -exponent c 1 . An increase in the stretching level has a similar effect in accentuating model nonlinearity (Fig. 9), but a careful tissue handling is required in order to avoid damage. The inclusion of a purely elastic term in the model is essential for modeling creep behavior, as seen from the inadequate shape of the Maxwell type of models. However, the fractional Maxwell model appears to be able to capture a representative behavior across all other tests. This indicates a possible redundancy of the elastic term across the other kinematics, which could lead to non-unique parameter identification. For example, a previous study involving the fractional Kelvin-Voigt model [23] indicated that appropriate torque response could be achieved with either a combination of large b 0 and α b values or small α b and b 0 = 0. Although this could be a consequence of employing similar forms for S e and S Q in both this study and [23], the inclusion of a purely elastic term is essential and could be well characterized through a combination of tests. In the above example, a relaxation test would have helped in identifying the b 0 and α b balance, as a small relaxation time (large α b ) would have led to a rapid decay towards the elastic limit b 0 S e , whereas a small α b would have shown a slow decay. Apart from the stress decay in relaxation tests, fractional parameters α a and α b also modulate the hysteresis levels, hence cyclic tests are required in determining their contribution. Generally it should be considered that α a ≤ α b , as otherwise the loss modulus shows non-physiological trends (Fig. 11).
We note that Bagley and Torvik [10] show, for a linear fractional Zener model, that α a = α b to satisfy the second law of thermodynamics under cyclic load. While non-dissipative behavior was not observed in the tests presented, for the case where (α a , α b ) = (0.2, 0.5), a 1 = 0.8, b 0 , b 1 , c 0 = 1 and frequencies below ≈0.0025 Hz, hysteresis amounted in a small accumulation of energy as opposed to a dissipation. Although this only occurs for small frequencies, this indeed suggest that α a = α b or α a = 0 is necessary to satisfiy the laws of thermodynamic.
The storage and loss moduli investigation over an extensive frequency range can also point towards the appropriate balanced contribution of a 1 and b 1 , as a 1 can modulate the frequency at which these moduli stop increasing according to a power-law, and b 1 can tailor the maximum level achieved. From this analysis it can be inferred that for a more complex model, like fractional Zener, an extensive testing is required in order to be able to separate parameter contribution -an oscillatory test at large stretch levels, creep, relaxation and frequency range investigation. If the latter is not available, then the use of the fractional Kelvin-Voigt model might be preferable, to avoid α a − α b ambiguity.

Limitations
In this study some assumptions and simplifications had to be adopted, due to the extensive amount of data produced by parameter variation. Firstly, the models were limited to be isotropic. Anisotropic forms would have shown a more pronounced differences between uniaxial and shear responses, whereas a redundancy of the two modes of deformation was observed here. Nonetheless, based on the knowledge gained in this study, a model with predefined set of feasible parameters (e.g., α a ≤ α b and a 1 ≤ b 1 ) can be extended to incorporate anisotropic forms for which the material parameters can be determined from further studies. Secondly, here the same underlying form was assumed for S e and S Q , with c e 1 = c v 1 . These forms could be altered, but the gain would not have been significant for the purpose of this study, whose focus is on the viscoelastic behavior exhibited by different model forms, rather than on the underlying hyperelastic behavior. Moreover, other types of material models could have been employed in lieu of the Fung-type model. For example, a study considering fractional viscoelastic models with Mooney-Rivlin, Fung and Ogden forms for a limited number of rheological tests is done in [23]. Nevertheless, considering multiple forms here would have further complicated the analysis of results, which currently consists of 4,216 models compared across 144 testing conditions. Model complexity could also be increased, e.g., by incorporating additional viscoelastic terms (b 2 D αc t S Q ). However, as the fractional Zener model showed enough versatility to capture typical tissue behavior across several testing protocols, additional terms might have simply led to redundancy in observed behavior. However, this idea might be pursued in combination with increasing material complexity, such as anisotropy.
The visual results presented in this study had to be narrowed down to several representative cases. More parameter combinations results are available in the Supplementary Materials, but an exhaustive representation would have been impossible due to the large number of curves produced. Nonetheless we believe that the discussion points addressed in Sects. 3 and 4 offer a good overview on models' capabilities and parameter influence.

Conclusions
Viscoelastic modeling of biological tissues is an area that has been explored to some extent in literature, but lags behind hyperelastic modeling amongst modelers. Classical viscoelastic models have been more readily employed, while fractional formulations have only been emerging more recently. The study presented here offers a systematic overview of nonlinear classical Maxwell, Kelvin-Voigt, Zener and their fractional model adaptions across a range of typical rheological tests. Different kinematic tests have been imposed (oscillatory uniaxial and shear deformations, constant load and constant strain), as well as variations in deformation level and frequency response. A comprehensive model parameter analysis was also performed across a broad spectrum, revealing the specific influence of each parameter on model behavior. While the results could not be presented exhaustively, due to the large number of curves produced, the relevant information was inferred by computing a similarity metric of each model-parameter combination against all the other combinations, or against a representative curve. Therefore, we could identify the scenarios which yield typical tissue responses, as well as understand the limitations of the other cases. As such, the findings of this study can be used in making an informed model selection when aiming to capture particular rheological traits. Alternatively, the insights gained here can help devising testing protocols such that model parameters can be well identified from varying tissue response.
Choosing an appropriate viscoelastic model depends on the behavior that needs to be replicated. If a simple, uniaxial cyclic stress response has to be characterized, then the classical Zener model would suffice. At low frequencies, a further simplification to classical Kelvin-Voigt model would work (i.e., a 1 = 0), whereas for intermediate to high frequencies, the simplification should be done towards the Maxwell form (i.e., b 0 = 0). Similarly, if relaxation and creep responses are not of interest, then one can select a simpler model like Maxwell or Kelvin-Voigt in order to replicate other types of behavior (e.g., cyclic). However, for realistic responses across a wide range of tests and frequencies, the fractional Zener model is the most robust, if employed with a reasonable set of parameters. This model could always be considered as a general starting point, with simplifications being carried as appropriate, depending on frequency range of interest and deformation type -cyclic or constant stress/strain. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.