Fractional viscoelastic models for power-law materials

Soft materials often exhibit a distinctive power-law viscoelastic response arising from broad distribution of time-scales present in their complex internal structure. A promising tool to accurately describe the rheological behaviour of soft materials is fractional calculus. However, its use in the scientific community remains limited due to the unusual notation and non-trivial properties of fractional operators. This review aims to provide a clear and accessible description of fractional viscoelastic models for a broad audience, and to demonstrate the ability of these models to deliver a unified approach for the characterisation of power-law materials. The use of a consistent framework for the analysis of rheological data would help classify the empirical behaviours of soft and biological materials, and better understand their response.


Introduction
Characterising and understanding how a material deforms when subjected to external forces is critical for many branches of industry, from food and material processing to product design. For example, in additive manufacturing, solid polymer filaments are melted and then extruded through a nozzle. The material flow behaviour during extrusion affects both the processing time required, and the strength of the final printed object [1,2]. In the food industry, the texture and mouthfeel of bread is largely dependent on its mechanical properties, which themselves are largely determined by the manufacturing process [3,4]. Furthermore, the mechanical response of polymers is often used to infer their composition or microstructure [5,6]. The scientific study of material deformation is known as rheology. To further its inquiry, rheology uses mathematical modelling to capture an approximate representation of the behaviour of materials. This facilitates the classification and comparison of materials, and enables predictions that can then inform engineering design choices and improve manufacturing processes.
Studying how living systems respond to mechanical stimuli is fundamental to gaining a comprehensive understanding of their functions and underlying cellular processes. It is increasingly recognized that the mechanical properties of cells have a great impact on human development and disease progression [7][8][9]. Alterations in the mechanical behaviour of soft tissues when subjected to external stimuli are often associated with diseases [10][11][12][13][14]. Being able to measure and quantify such changes can provide insights into disease evolution, which in turn guide advancements in diagnostic tools and treatments [15][16][17]. In tissue engineering, knowing the response of human tissue to external forces is key to developing replacements that restore the structure and functionality of damaged tissue by matching the properties of the surrounding environment, which promotes cell in-growth and reduces the risk of implant failure [18][19][20][21].
There are two common limits in the behaviour of materials. Elastic materials are solids whose mechanical stress is related to their deformation state with respect to a reference, stress free, geometry; the work done to deform them is entirely stored as recoverable elastic energy. Newtonian Semi-hard cheese Pure ice Asphalt sealants Epithelial cell Figure 1: Stress relaxation response of several materials follows a power-law behaviour. Examples of power-law materials include: (a) common polysaccharide (Xanthan gum) used as food additive [22], bread dough [23], synthetic polymers such as nylon of diameter 1.125mm [24], single collagen fibrils [25]; (b) semi-hard zero-fat cheese [26], pure ice at -35 • C [24], asphalt sealants [27] and single MDCK epithelial cell [28]. For clarity the data has been separated into two subplots and the stresses have been normalized by the initial value for each data set.
fluids are viscous materials whose rate of deformation is function of their mechanical stress; the mechanical work required to deform the material is fully dissipated. Most materials however exist on a spectrum between the elastic and Newtonian limit behaviours; they are referred to as viscoelastic materials. A defining characteristic of viscoelastic materials is their time-dependent behaviour. When a constant deformation is applied to viscoelastic materials, their internal stress decreases with time -this process is known as relaxation. When viscoelastic materials are subjected to a constant load (stress), their deformation (strain) increases with time -this process is known as creep. In applications where long-term durability and shape stability are required, creep may not be a desirable material property. In contrast, creep may be a desired property in manufacturing processes where a material is extruded or flows across cavities -e.g. in additive manufacturing. Various empirical methods are available to quantify the response of viscoelastic materials. Relaxation tests (in which a constant strain is applied and stress is recorded) and creep tests (in which a constant stress is applied whilst strain is recorded) constitute two testing paradigms for viscoelastic materials. Another common testing methodology requires the application of an oscillatory load to the material. Typical instruments to collect such data are Dynamic Mechanical Analysers (DMA) for materials with solid-like behaviours at long time-scales and rheometers for liquid-like materials that need to be studied in shear deformations only. Relaxation, creep and oscillatory tests are the three canonical viscoelastic testing paradigms; the most appropriate test can depend on several factors including the experimental hardware available, and the most relevant physical modes of deformation for the application or research question. Mathematical models are commonly used to describe the creep, relaxation and oscillatory behaviour of viscoelastic materials. One of the primary aims of modelling is to extract model parameters and relate them to the underlying molecular or microstructural deformation mechanisms. In other contexts, such as the food industry, the parameters can be related to sensory perception. Nevertheless, the identification of a suitable model facilitates comparison between the behaviour of different materials, states of organisation or environmental conditions. The theory of linear viscoelasticity is a widely used approach to analyse experimental data that yield to a well-defined mathematical representation of the stress-strain-time relation -often referred to as constitutive relationship. By employing those mathematical methods it is possible to extract parameters that uniquely describe a given material -often referred to as material properties -and allow the prediction of the evolution of stress and strain within the material under arbitrary loading scenarios.
Advances in mechanical testing hardware have revealed that many materials possessing a com-plex microstructure exhibit a characteristic power-law signature in their creep, relaxation and spectral behaviours. A selection of data from the literature is shown in figure 1. Further examples include biological materials -e.g. tissues [29][30][31], single cells [28,29,[32][33][34][35][36], intra and extra cellular components [25,37,38]; gels [39,40], polymers [24,41], concrete [42,43], asphalt [44][45][46], ice [24], and food-e.g. cheese [26,47], dough [23,48]. This behaviour represents a challenge to the commonly used viscoelastic models in obtaining a unique mathematical description of the constitutive relation with a limited number of parameters. The power-law behaviour can be approximated with traditional viscoelastic models via a large number of model parameters, which greatly hinders physical interpretation. To circumvent this limitation, traditional models are often discarded in favor of empirical fitting functions: mathematical ansatzes derived from qualitative inspection of the data rather than a governing constitutive equation. This approach precludes comparison of parameters across research studies, which can lead to a multitude of different interpretations of similar phenomena [49]. A promising solution relies on the use of fractional calculus, a branch of mathematics that extends integration and differentiation operators to non-integer order [50,51], to enrich the classical linear viscoelastic framework [52,53]. This led to the development of a new formalism for the modelling of viscoelastic materials known as fractional viscoelasticity. Fractional viscoelasticity has been applied to complex geological and construction materials such as bitumen (asphalt) [44,45], concrete [42,43], rock mass [54][55][56][57][58][59], waxy crude oil [60,61], as well as polymers and gels [41,[62][63][64][65], and food [26]. Numerous examples can also be found of fractional viscoelasticity applied to biological materials such as epithelial cells [66], breast tissue cells [67,68], lung parenchyma [69], blood flow [70,71], as well as red blood cell membranes [72]. In spite of the examples reported above, when we consider the number of power-law responses observed in the literature (some of which are modelled empirically), fractional viscoelasticity remains significantly underused. By using a consistent formalism for the analysis of power-law materials, a direct comparison of parameters across studies can be achieved [66], which greatly extends the use of available data.
The main aims of this review are to demonstrate the applicability of fractional viscoelasticity to a wide range of real problems, and to carefully illustrate its advantages over commonly used modelling approaches. Furthermore, we aim to present fractional viscoelasticity in a way that is accessible to researchers without specialised knowledge of fractional calculus. To achieve this, we first illustrate the limitations of spring-dashpot models in capturing power-law viscoelasticity. Following this, the mathematical foundations of fractional viscoelasticity are introduced. We will describe the behaviour of networks of fractional viscoelastic elements with a focus on their limit behaviours are short and long time-scales. A number of studies which originally used empirical or spring-dashpot based traditional models are then revisited in order to demonstrate the benefits of fractional viscoelastic models. A physical interpretation of the phenomena is discussed where possible.

Introduction to linear viscoelasticity
The linear theory of viscoelasticity provides a powerful mathematical framework to link stress, strain and time and provide predictions of stress and strain distribution during arbitrary loading conditions. A material is assumed linear viscoelastic if the stress function σ(t) and strain function (t) are linearly connected: (i) if the strain function (t) is multiplied by a constant factor, the resulting stress would be scaled by the same constant, and vice-versa. (ii) the response to a strain (or stress) that is a linear combination of two arbitrary strain (or stress) functions is given by the same linear combination of their two individual responses. Most materials exhibit linear, or quasi linear, behaviour for small deformations while their response becomes nonlinear at large deformations. The linear theory of viscoelasticity is a commonly used approximation for the study of materials' behaviour as it leads to an accessible and manageable mathematical formulation of the stress, strain and time relationship. Furthermore, the theory of linear viscoelasticity represents a good starting point for the study of more complex nonlinear responses that often gives rise to more complex models [73][74][75][76].
A common test performed to characterise a viscoelastic material is to analyse the stress response of the material when subjected to a constant strain (relaxation test). The implication of the linearity assumption is that the resultant stress function scales linearly with the magnitude of the step in strain [77]. Thus, for a step in strain of amplitude 0 at time t = 0, the resulting stress σ(t) can be written where G(t) is the relaxation modulus, a monotonically decreasing function. Similarly, in the case of a creep test, the strain response (t) is proportional to the stress step amplitude σ 0 that is imposed at t = 0: where J(t) is the creep modulus, a monotonically increasing function of time.
A key consequence of the linearity assumption is that the response of the material at time t is given by the sum of the responses to the perturbations imposed at all previous times. Therefore, given an arbitrary stress or strain history, we can compute either the strain response or stress response respectively by integrating the entire history of infinitesimal 'step' excitations [78,79]. Assuming that the material does not age, i.e. that its mechanical behaviour does not change with time during the course of a measurement, the strain (or stress) response at time t of a step in stress (or strain) at time τ would be given by the relaxation function G(t − τ ) (or creep function J(t − τ )). This results in the following convolution integrals: Another consequence of linearity is that, when a viscoelastic material is subjected to a sinusoidal stress, the recorded strain has the same frequency, but with a phase difference that may depend on frequency. For purely elastic materials the phase difference is 0°and for purely viscous materials the phase difference is 90°. For viscoelastic materials, the phase difference lies between these two values. If we assume a linear material, by considering an oscillatory excitation (t) = e iωt , we obtain a stress response σ(t) = G * e iωt [79] from which we can define the so-called complex modulus, or dynamic modulus: where ω is the frequency. The real part of this complex stress response, G (ω), is defined as the storage modulus (as the energy is stored in an ideal elastic material). The imaginary part of the response, G , is defined as the loss modulus (as energy is dissipated in a purely viscous material). For a linear material, the three moduli, relaxation, creep and dynamic moduli, are directly related. Their relationship can be expressed in the Laplace domain. For instance, G(s) J(s) = s −2 , where G(s) and J(s) are the Laplace transforms of equations (3) and (4) respectively [73]. Therefore, we can predict one type behavioural mode (e.g. creep) from observing a different behavioural mode (e.g. relaxation).

Characteristics of traditional linear viscoelastic models
The expressions above are generic and applicable to any material that behaves in a linear manner, being solid and liquid like. Specific expressions for the moduli -i.e. relaxation, creep and dynamic moduli -can be deduced in several ways. A common approach relies on the use of fundamental viscoelastic units that can be combined in series or parallel -this element-by-element model building process is analogous to electrical circuit models. The most commonly used viscoelastic units are the Hookean spring and the Newtonian dashpot, with governing equations of σ(t) = k (t) and σ(t) = η˙ (t) respectively. The simplest models that can be formed from these two viscoelastic units are their series (Maxwell model) and parallel (Kelvin-Voigt model) combinations, both of which are shown in figure 2. The creep and relaxation moduli can be derived from the corresponding differential equation by: a) assuming that the stress is constant and solving for the differential equation in strain, leading to the creep modulus, and b) assuming that the strain is constant and solving for the differential equation in stress, leading to the relaxation modulus. The complex modulus is derived from the constitutive differential equation, either by considering a oscillatory forms of the stress and strain functions, or by taking the Fourier transform of the differential equation, leading to an algebraic equation which can be rearranged for G * =σ(ω)/ˆ (ω), where the hat (ˆ) symbol denotes a Fourier transformed function. As the differential equations relating stress and strain of traditional spring-dashpot models take the form of a linear ordinary differential equation with constant coefficients, both relaxation and creep moduli involve exponential functions (see figure 2). As an example of this, consider the stress response of a Maxwell model to a step in strain: where G(t) is the relaxation modulus. The ratio τ = η/k has the units of time and can be considered as the characteristic relaxation time of the material. This time-scale is a useful quantity; it represents the time required for the stress to fall to 1/e of its initial value. More importantly, the time-scale may describe a physical relaxation process more effectively than the values of η and k separately, such as rates of binding or unbinding of molecular structures. To capture more complex material behaviours, additional elastic springs and viscous dashpot elements may be used. These additional elements give rise to multiple time-scales. Some models, for example the generalised Maxwell model shown in figure 3 (a), give rise to an arbitrarily large number of time-scales. (See [73,80,81] for more details about spring-dashpot based linear viscoelastic models.) Power-law viscoelastic materials correspond to systems with a broad range of relaxation/creep time-scales arising from dissipation/deformation modes occurring at different time and length scales. Mathematically, power-law behaviour consists not of a discrete number of time-scales, but a continuous distribution of time-scales. For this reason, power-law rheology can only be approximated by the spring-dashpot models that yield exponential terms with discrete time-scales [82][83][84]. The larger the number of exponential terms, the closer the approximation [85][86][87][88]. We demonstrate this in figure 3 by incrementally increasing the number of Maxwell arms in a Generalised Maxwell model, each of which contributes a material time-scale. With four Maxwell arms, a good approximation of a power-law is obtained over the time domain of the data (figure 3 (b)). Another informative perspective on the approximation process can be gained via relaxation spectrum, which represents the relative behavioural dominance at different times. The relaxation spectrum found here for the power-law relaxation modulus was approximated as [79,89] where G(t) is the relaxation modulus. Figure 3 (c) shows the relaxation spectrum of the powerlaw material in logarithmic scale, which is a power-law distribution of time-scales with the same exponent as the original relaxation data. The optimal fit of the various tested Maxwell models can be seen to yield time-scales that are equally spaced in logarithmic scale. As more Maxwell arms are added, a broader coverage of the true power-law distribution is attained. This qualitative assessment of the above relaxation responses and their corresponding spectra shows that spring-dashpot models might be sufficient to approximate power-law behaviour. However, the approach has a number of disadvantages. A large number of model parameters make the computational analysis far more expensive, and interpretation of their physical meaning becomes significantly more difficult. The characteristic times extracted from the fitting process mostly capture the time window of the data and do not represent intrinsic material properties. Moreover, the fitting would optimise the match with data in the region fitted, but provide poor predictions at time-scales shorter of longer than those provide in the original data.
To circumvent the above disadvantages inherent to spring-dashpot based approximations, empirical power-law expressions have been used as a simple way to capture experimental measurements [29,32,33,90,91]. These expressions are essentially a mathematical ansatz, formulated to capture the main qualitative features of the experimental data. The end result of this formulation is usually a single modulus of interest that may or may not be analytically converted into the other moduli, necessitating numerical analysis for the prediction of the material response to other external loading patterns. Further, the use of such ad-hoc models limits the scope of the measurements since the model parameters may not be easily comparable across studies.

Fractional viscoelastic models
With regards to the modelling of power-law viscoelastic materials, we have now outlined the advantages and disadvantages of the phenomenological spring-dashpot approach and the empirical ad-hoc approach. In this section we introduce an additional viscoelastic element, the springpot, which is able to capture a power-law behaviour with a minimum number of parameters. We then demonstrate how the springpot can be combined with other elements in series and parallel configurations to capture complex power-law signatures in various contexts.

The springpot captures power-law behaviour
By inspecting the power-law relaxation data (figure 1), it can be seen that the relaxation modulus takes the form G(t) = At −β , where A is a constant and β the power exponent. Substituting this into equation 3, the relationship between stress and strain for a power-law material becomes: A branch of Mathematics called Fractional Calculus provides the definition of the generalization of the differentiation operation to non-integer order as [92] where Γ(·) is the gamma function. Redefining the coefficient A = c β /Γ(1 − β) in equation 8 yields the following result: which provides a simple constitutive equation for power-law materials. A more detailed derivation of this relationship is presented in the next section. The springpot's schematic symbol and relationships to the spring and dashpot elements are shown in figure 4. For dimensional consistency, the unit of the constant c β must be Pa s −β . Due to its unusual physical unit, the parameter c β lacks a tangible physical meaning, but is often interpreted as the 'firmness' of a material [93]. It should be noted that in the limiting cases where β = 0 or β = 1, the springpot reduces to a spring or dashpot respectively, and the constant c β represents the elastic spring constant, k (Pa), or the dashpot viscosity, η (Pa s), respectively.
To overcome the use of a constant c β whose physical meaning may be unclear, the constitutive equation may be written as σ(t) = λ 0 τ β 0 d β (t)/dt β , where λ 0 has units of Pa and τ 0 is a characteristic time that has units of s [94,95]. However, it is not possible to design experiments to separately measure λ 0 and τ 0 . Therefore, in practical problems, it is more relevant to define the fractional element with two parameters c β and β [96], and this is what we use in the current review.

Derivation of the springpot's governing equation
Although not essential for springpot practitioners, it is illuminating to derive the springpot's governing equation. Whilst the correspondence between power-law viscoelasticity and the fractional derivative has been presented before [97,98], here we attempt summarize the main steps of the derivation in a pedagogically accessible way.
Fractional derivatives can be understood as generalised integrals. Integrating a function (here (t)) an integer n number of times, we get: (11) where the final simplification on the RHS made use of Cauchy's repeated integral formula. Equation 11 can in fact be generalised to any positive real number α by recalling the extension of factorial numbers to non-integer values via the Gamma function Γ(α) = (α − 1)!
If we substitute α = 1 − β into equation (12) and f (t) = d (τ ) dτ , we can then re-write the expression for the total stress reported in equation (8) as follows Since the fundamental relation (I a+b f )(t) = I a I b f (t) holds [92], equation (13) can be expressed as where I −β (t) is the definition of fractional derivative D β (t). Equivalent to the above, we can write which is the governing equation of the springpot [53].
It should be noted that in the literature, different definitions of the fractional derivative have been put forward. Here we use the Caputo's derivative definition since it naturally arises from experimental observations and it has shown to have a better applicability to real problems, where initial conditions are known in terms of derivatives of integer index [99]. If we assume that the system is at rest for time t ≤ 0, the fractional derivative is given by [100] where 0 < β < 1. So to obtain the value of the Caputo fractional derivative at time t, we must integrate from the initial time t = 0. In other words the fractional derivative is a non-local operator. This is in contrast to integer order derivatives, whose value is determined by the limiting behaviour of the function at the evaluated point. Although hysteresis is also present in the generic hereditary integral shown in equation 3, it is a fundamental feature of the fractional derivative and consequently the springpot. Lastly, it is important to note that the fractional operator is a linear operator. Therefore, the springpot element lies within the framework of linear viscoelasticity. equation (15) and making use of the definition of the fractional derivative in equation (16), we can extract the relaxation modulus G(t) of the springpot [99]: The relaxation modulus is as expected a power-law function of the time whose exponent matches the order of the fractional derivative. As a result, the corresponding relaxation spectrum of the springpot is also a simple power-law of exponent β, and therefore captures very well, by design, the behaviour of power-law materials, with only two parameters. The creep modulus J(t) can be deduced from the expression of G(t). Both are related in the Laplace domain by G(s) J(s) = s −2 .
The Laplace transform of G(t) in equation (17) is given by G(s) = c β s β−1 . Therefore, J(s) = 1/(c β s 1+β ), whose inverse Laplace gives the creep modulus [99] Qualitatively, the linear elastic solid and the Newtonian fluid behaviours are the two limit behaviours of the springpot for β = 0 and 1 respectively. To understand the interpolation of the springpot between elastic and viscous behaviour, we plot the springpot response to a step in strain for different values of β ( figure 5 (a)). When a dashpot (β = 1) is subjected to a step in strain, the stress is initially infinite but immediately dissipated afterwards ( figure 5 (a), blue curve). By gradually decreasing the exponent β the time taken by the springpot to dissipate the stress increases until the limit of the linear elastic solid is reached. When a spring (β = 0) is subjected to a constant deformation, the stress reached is proportional to the spring constant k ( figure 5 (a), red curve).
To illustrate the creep behaviour of the springpot, we consider its response to a step in stress σ 0 imposed for a limited time t * -loading phase -and then returned to 0 -unloading phase. Figure 5 (b) shows plots of the system's response for different values of β during the loading and unloading phases. A spring generates a strain directly proportional to the stress, as expected, and it immediately returns to its original length upon unloading. A dashpot linearly deforms while the stress is imposed, and does not return to its initial state after unloading. The response of the springpot is more complex. By applying the superposition principle, we can write an expression    for the springpot response (0 < β < 1) after the load is removed at time t = t * that is given by noting that when t tends to infinity, the strain tends to zero. Therefore, when the stress is removed, the strain in a springpot is always completely recovered; the springpot has shape memory, although energy would be dissipated during deformation. The greater the power-law exponent of the springpot β, the larger the time required for total strain recovery, as shown in figure 5 (b). For example, by looking at the response of a springpot with β = 0.7 (green curve in figure 5 (b)) over a finite time interval, one might qualitatively interpret the slow recovery process as evidence for plastic behaviour -although this is erroneous because strain eventually converges to zero. This shows that long experiments are needed to fully capture the rheology of power-law materials. By taking the Fourier transform of equation (15), the complex modulus of a springpot can be found where ω is the frequency. Separating the real and imaginary parts, we find the storage and loss moduli respectively They follow a power-law behaviour with the same exponent β. When β = 0 the loss modulus is zero as in the case of the spring, whilst when β = 1 the storage modulus is zero as in the case of a dashpot ( figure 6 (a)). The storage and loss modulus exactly match when β = 0.5 ( figure 6 (a)). For β < 0.5 the storage modulus is always greater than the loss modulus, whilst the opposite is true for β > 0.5 ( figure 6 (b)). The phase angle δ between the excitation and the response is related to the storage and loss moduli by tan(δ) = G /G , from which we can derive the retardation phase for a springpot as δ = π 2 β, constant for all frequencies.

Examples of the use of the springpot in practical cases
The first material investigated using fractional calculus was bitumen in 1944, in a study conducted by Scott-Blair and Veinoglou [44]. Since then, use of the springpot has been somewhat erratic, possibly because of its non-trivial mathematical foundations. However, the springpot has found use in a diverse array of materials outside of the geological and construction-materials contexts; for example, in gels. An early example is the work of Winter and Chambon [84] who derived a springpot-like modulus for crosslinking polymers at their gelation point, which was used for analysis of polydimethylsiloxane gel data. More recently, a springpot has been used in modelling asphalt mixtures [45] and to capture the oscillatory rheological response of agarose, a polysaccharide extracted from red algae [62]. More examples of springpot usage can be found in the context of biological materials more generally. One of the first applications of the springpot in tissue biomechanics was for the study of the viscoelastic properties of lung [101]. More recently, it has been used to capture the relaxation behaviour of human arteries [102] and the dynamic (oscillatory) response of brain tissue for the understanding of neurological disorders [103]. There are also several biomechanical studies which make use of a power-law empirical function that corresponds to one of the three springpot moduli discussed above. For example, the empirical function used to analyze the creep response of single cells can be easily related to the springpot's creep modulus [32,104,105]. Similarly, an empirical function used to analyze the relaxation response of smooth muscle cells is equivalent to that of a single springpot [106].

Generalised fractional viscoelastic models
Many materials exhibit power-law viscoelastic behaviour. However, the power-law regime is often limited in time; some materials show a power-law response at short time-scale that converges to a well-defined plateau (solid like behaviour) [29,72], whilst others may exhibit a power-law response followed by a continuous flow of the material (fluid like behaviour) [107,108]. In order to capture the diversity of behaviours, ad-hoc empirical moduli and large networks of springs and dashpots (i.e. exponential terms) have been often used. However, as previously discussed, these approaches each have disadvantages.
A third useful approach is to combine springpots with themselves and with traditional spring and dashpot elements. In this way, the many advantages of the spring-dashpot network approach can be extended to concisely capture power-law viscoelastic regimes. Such models are often referred to as generalised viscoelastic models and they have been mathematically characterized in the past [94,95,109]. The two simplest configurations are two springpots in series or in parallel, which are the fractional analogues of the Maxwell and Kelvin-Voigt models respectively. Examples of the use of the generalized models in practical cases are reported in table 1 (see [110] for a further review of fractional viscoelasticity in geotechnical engineering). In the following section, we demonstrate the behaviour diversity that these generalised models are capable of.

The Fractional Maxwell and Kelvin-Voigt Models: Heuristic Overview
The choice of which viscoelastic model to use is largely informed by the qualitative behaviour of the experimental data. A key advantage of schematically representing viscoelastic models as networks of elements is that it facilitates visual intuition for a model's behaviour. Despite the mathematical complexity of the springpot, limit behaviours can be obtained as commonly done for springs and dashpots models. To address this, we provide below insight on the qualitative behaviour of the two simplest generalised fractional viscoelastic models, which can be extrapolated to more complex models.
When two springpots of exponents α and β (with α > β) are placed in series, the element with the lower exponent (more "springy") dominates the short time-scale response, while the springpot with the higher exponent determines the long time-scale behaviour (see figure 7 first column). This is true for both the relaxation and creep responses, thus When two springpots are combined in parallel the opposite is true: both relaxation and creep responses are dominated by the springpot with higher exponent at short time scale, while the springpot with the lower exponent dominates at long time-scale (see figure 7 second column), thus Indeed, traditional viscoelastic models can be seen as special cases of the generalised fractional models obtained by specialising the springpots to either springs or dashpots. The qualitative behaviours reported above are consistent with those observed in conventional viscoelastic models, where the Maxwell model shows a liquid-like behaviour at long time scale (dominated by the dashpot, α = 1), while the Kelvin-Voigt model shows a solid-like behaviour at long time-scale (dominated by the spring, β = 0). Although the above heuristics directly apply to the two-springpot models, the insights can be extrapolated to more complex models in a straightforward manner. For example, consider the Standard Linear Solid model (figure 2) but with the spring k 1 replaced by a springpot. If we look at the global network, the k 2 spring in parallel dominates the long-time scale response given its lower power. However, if we focus only on the upper arm, the model will follow the series configuration behaviour. Thus, the short time-scale response is dominated by the springpot and the intermediate response by the dashpot. A more graphical description of this heuristic process has been recently presented by Bonfanti et al. [66]. To build further intuition, readers may also refer to the Annex, where the relaxation and creep behaviours of more complicated models are Table 1: Examples of practical uses of generalized models.

Fractional Kelvin-Voigt
Modelling of human prostate tissue to develop novel criteria for cancer detection [111].
Modelling of the oscillatory response of canine liver [112].
Modelling of the relaxation response of beast cells and tissue samples [40].
Modelling of the tissue mimicking materials CF-11 and gelatin [39].
Modelling of the viscoelastic response of potato starch gel [113].
Modelling of a three-dimensional particulate gel [63].
Studying of the applicability of the model to wellbore creep [114] Fractional Maxwell Modelling of the creep behaviour of rocks [55].
Modelling of the viscoelastic response of tight sandstone [56].
Modelling the properties of food gels [26].
Modelling the mechanical properties of collagen gel [117].
Fractional Standard Linear Solid model Study of single red blood cells [72].
Modelling of breast cancer cells to develop new diagnostic tool [68].
Modelling of artery walls for the study of aneurysms [70,71]. Quantification of changes in viscoelastic response of lung parenchyma due to trauma [69].
Modelling of the oscillatory response of cancerous cells; the information was then used to selectively attack malignant cells during treatments [118].
Modelling of the viscoelastic response of brain tissue [119].
Modelling of the viscoelastic properties of the glassy amorphous polymer polymethyl-methacrylate [41].
briefly reported.
In the following sections we demonstrate the benefits of using fractional models in a wide range of practical applications. We first show that a number of empirical functions previously introduced in the literature actually possess an equivalent fractional model. We then demonstrate how the use of fractional models often allows us to capture the mechanical behaviour of power-law materials more accurately than spring-dashpot viscoelastic models while using less parameters.

Examples of empirical functions equivalent to fractional viscoelastic models
Structural damping model. A number of dynamic (oscillatory) viscoelastic tests on biological tissue were fitted with an empirically derived model known as the structural damping (or hysteretic damping) model [33,[121][122][123][124][125][126] whose complex modulus given by: where η = tan β π 2 is commonly called the hysteresivity or structural damping coefficient, β is the power-law exponent, G 0 and ω 0 are two scaling factors for stiffness and frequency respectively, and µ the Newtonian viscosity. This modulus can be shown to be exactly equivalent to that of a semi-fractional Kelvin-Voigt model consisting of a springpot in parallel with a dashpot. In fact, equation (24) can be re-written as: which is equivalent to the complex modulus of a dashpot in parallel to a springpot G * (ω) = c α (iω) α + c β (iω) β , when α = 1, c α = µ and c β = G 0 /ω β 0 .
Two power-law regime. The frequency response of single cells often exhibits two power-law regimes [38,120]. A lower-power exponent at low frequencies, followed by a higher power-law exponent at high frequencies. Such behaviour has been frequently been analyzed by the empirically derived superposition of two power-laws where one exponent is often fixed to β = 3/4. The expression above is exactly equivalent to the complex modulus of the Fractional Kelvin-Voigt model shown in the inset of figure 8 and consisting of two springpots in parallel. Examples of fitting the storage and loss moduli of single cells to the Fractional Kelvin-Voigt model are shown in figure 8.
Power-law cut-off. Recently, Khalilgharibi et.al. [29] studied the relaxation response of epithelial monolayers, a sheet of cells devoid of substrate, to uncover the subcellular components responsible for stress dissipation. The relaxation response consists of an initial power-law phase in the first 5 s, followed by an exponential phase that reaches a plateau at ∼ 60 s. This biphasic relaxation response was captured by extending the power-law empirical function to include an exponential regime More recently however, traditional viscoelastic elements have been combined with a springpot to characterize the viscoelastic response of epithelial monolayers at both short and long time-scale [66]. The novel configuration features a springpot that sustains all the load at short time-scale, after which the load is then slowly transferred to a dashpot placed in series which gives rise to the exponential behaviour. This dissipative branch is placed in parallel with a spring to obtain the final plateau (inset in figure 9 (a)). The mathematical expression for the relaxation modulus of the fractional model introduced by Bonfanti et al. [66] is where c β and β are the springpot parameters, η the dashpot viscosity, k the spring stiffness, and E a,b (z) is the Mittag-Leffler function, a special function that often arises from the solution of fractional differential equations [128,129]. By fitting the fractional model to the empirical functions we can demonstrate that they are approximately equivalent (figure 9 (a)) and therefore the empirical function can be mapped into the fractional viscoelastic framework. Interestingly, the same empirical model was independently developed to capture the relaxation response of one of the most widely used elastomers, polydimethylsiloxane (PDMS) [127], and the rheological behaviour of alginate-based gels [130,131]. The study of the liquid-solid transition during gelation has attracted significant attention in the past [132]. Friedrich et al. [127] reported the rheological response of stoichiometrically balanced polydimethylsiloxane (PDMS) in which the cross-linking reaction was interrupted at different times, before and after the critical point (referred to as the gel point). As shown in figure 10, the four-parameter fractional solid model can also be applied here to accurately captures the behaviour of the time-evolving cross-linking reaction preand post-gelation.
Physically, during gelation a polymer network forms. The critical gel point is associated with an abrupt change in viscosity [133], which is defined as the creation of the first percolation cluster that spans the sample. By observing the values of the fitted parameters reported in table 7 in Appendix, the fluid-like behaviour of the PDMS pre-gel is confirmed by the zero value of the stiffness of the spring in parallel. As expected from physical considerations, at the gel point we observe a rapid increase of the viscosity (figure 12 in Appendix). Post-gel the PDMS behaves as a solid, which is confirmed by the rapid increase of the stiffness k. As discussed above, the storage and loss moduli of the springpot are exactly equal when β = 0.5, or physically, when the springpot is exactly intermediate between a liquid and a solid. The physical relevance of this parameter is confirmed from the fitted springpot coefficient shown in table 7, where β at time t − t c = 0 (gel point) is ∼ 0.5.
Modified power-law models. Another empirical model, recently used for the analysis of the relaxation response of various benign and malignant cell lines, is the modified power-law (MPL) model [87] which can be written as where E 0 is the instantaneous or 'glassy' modulus, E ∞ is the plateau or 'rubbery' modulus, τ is a time scaling, and α determines the power-law gradient of relaxation. In contrast to a regular power-law model, the MPL model has the convenient property of being well defined at short time scales. Its behaviour is notably similar to a Fractional Standard Linear Solid model (which is equivalent to a Fractional Zener model with two springpots specialised to springs), as shown in figure 9 (b). The relaxation modulus of the Fractional Standard Linear Solid model is defined as where k β and k γ are the two spring constants, c α and α are the springpot parameters (with a prime added to α to avoid confusion with the MPL α parameter), and E is the Mittag-Leffler function [128,129]. The similarity between the two models was also discussed by Bagley [134], who made comparisons between their relaxation spectra. Here we show their similarity by directly matching their boundary condition parameters at t = 0 and t → ∞ and fitting the Fractional Standard Linear Solid parameters c α and α to the MPL ( figure 9 (b)). The fractional model fits the MPL well, especially at longer time scales where both the Mittag-Leffler function and the MPL asymptotically approach a simple power-law [135]. Interestingly, the MPL model has also been used in several asphalt and asphalt-concrete studies [46,136,137] indicating that the fractional standard linear solid may have utility in this field.

Fractional models capture complex power-law behaviours with less parameters
As discussed, spring-dashpot viscoelastic models can be used to capture power-law viscoelasticity but this results in an abundance of model parameters. We have shown that the springpot provides an efficient way to describe the rheology of power-law materials with only two parameters. To further illustrate the ability of fractional models to accurately capture complex power-law responses with a minimum number of model parameters, we now re-analyse the time-response of different power-law materials, originally modelled using spring-dashpot networks, by using fractional viscoelastic models. We first analyse the rheological behaviour of single cells presented by Darling et al. [34]. They investigated the relaxation response of zonal articular chondrocytes by the use of AFM (Atomic Force Microscopy). The Standard Linear Solid model originally used in the paper struggles to capture the short-time scale response (see figure 11 (a)). We propose the use of the fractional Kelvin-Voigt model (springpot in parallel to a spring) that involves the same number of parameters (see inset in figure 11 (a)). A cursory look at the comparison between the fits of traditional and fractional viscoelastic models to the relaxation response of zonal articular chondrocytes shown in figure 11 (a) reveals the ability of the fractional model to significantly improve the quality of the fitting by accurately capturing the fast relaxation response while using the same number of parameters as the spring-dashpot model used in the original paper.
Recently, Zhang et al. [60] tested tomato mesocarp cells under high speed microcompression to understand how industrial handling may affect the integrity of fresh fruits. To successfully capture the distribution of time-scales that gives rise to the power-law regime at short time-scale, a 2 time-scale Standard Linear Solid that requires five parameters was originally used. By using a fractional Maxwell model (springpot in series to a dashpot) we achieve the same excellent fit to the data whilst reducing the number of parameters to three ( figure 11 (b)).
The two-time-scale Standard Linear Solid model has recently been applied to the analysis of the rheological behaviour of other materials, such as polycaprolactone bioactive glass tested under compression [139] and single collagen fibrils from the extracellular matrix under tensile testing [25]. Given that they present the same qualitative behaviour as the epithelial monolayers (power-law followed by exponential behaviour until a final steady-state value is reached), we have successfully applied the fractional model recently developed by Bonfanti et al. [66] to data from these other materials. The model accurately fits both short and long time-scale responses while using one less parameters (figure 11 (c)-(d)).

Conclusions
This review has demonstrated the ability of fractional viscoelastic models to accurately capture the rheological responses of a broad range of materials, identifying materials parameters that account for the wide distribution of time-scales involved in power-law behaviour. Despite the limitation to the linear response, fractional models exhibit rich behaviours consistent with empirical data in both relaxation and creep experiment. Looking at large deformations and failure of complex materials remains however a challenge that cannot be tackled by linear rheological models and we anticipate significant efforts to tackle these questions in future.
The main impediments to the use and dissemination of fractional viscoelastic models appears to be their relatively intricate mathematical formalism and the difficulty of their numerical implementation. This document illustrated the qualitative response of the springpot element and its composition into simple networks in series and in parallel configurations, helping the reader to build intuition. An annex to the review [REF HERE] provides an exhaustive list of rheological models including up to three elements, alongside the analytical expression of their moduli and graphs illustrating their behaviour. We also developed a software library to facilitate the fitting of experimental data and prediction of power-law behaviours [140]. The software package allows non-experts to fit their own data using a wide selection of viscoelastic models, accounting for complex loading patterns (all figures of this work have been created using RHEOS). Altogether, these elements significantly lower the barriers to using fractional models to analyse the rheology of power-law materials.
Although fractional models extend the range of behaviours that can be modelled with rheological elements, they do not provide on their own explanations for the underlying mechanisms giving rise to the observed macroscopic power-law behaviour. A deeper understanding of soft material mechanics would require a more systematic theoretical analysis of experimental data that would provide a physical underpinning for the emergence of power-law behaviours. However, fractional models provide systematic approaches to capture material parameters that can be compared across studies, and such comparison is likely to provide a better handle on the underlying physical significance of power-law behaviours.