DOI:
10.1039/D0SM00354A
(Review Article)
Soft Matter, 2020,
16, 60026020
Fractional viscoelastic models for powerlaw materials†
Received
28th February 2020
, Accepted 30th May 2020
First published on 8th June 2020
Soft materials often exhibit a distinctive powerlaw viscoelastic response arising from broad distribution of timescales 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 nontrivial 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 powerlaw 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.
1 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.^{14–16} Alterations in the mechanical behaviour of soft tissues when subjected to external stimuli are often associated with diseases.^{17–21} Being able to measure and quantify such changes can provide insights into disease evolution, which in turn guide advancements in diagnostic tools and treatments.^{22–24} 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 ingrowth and reduces the risk of implant failure.^{25–28}
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 fluids are viscous materials whose rate of deformation is a 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 timedependent 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 longterm 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 solidlike behaviours at long timescales and rheometers for liquidlike 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 welldefined mathematical representation of the stress–strain–time relation – often referred to as a 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 complex microstructure exhibit a characteristic powerlaw signature in their creep, relaxation and spectral behaviours. A selection of data from the literature is shown in Fig. 1. Further examples include biological materials – e.g. tissues,^{29–31} single cells,^{13,29,32–36} intra and extra cellular components;^{10,37,38} gels,^{39,40} polymers,^{9,41} concrete,^{42,43} asphalt,^{44–46} ice,^{9} and food – e.g. cheese,^{11,47} dough.^{8,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 powerlaw 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}

 Fig. 1 Stress relaxation response of several materials follows a powerlaw behaviour. Examples of powerlaw materials include: (a) common polysaccharide (xanthan gum) used as food additive,^{7} bread dough,^{8} synthetic polymers such as nylon of diameter 1.125 mm,^{9} single collagen fibrils;^{10} (b) semihard zerofat cheese,^{11} pure ice at −35 °C,^{9} asphalt sealants^{12} and single MDCK epithelial cell.^{13} For clarity the data has been separated into two subplots and the stresses have been normalized by the initial value for each data set.  
A promising solution relies on the use of fractional calculus, a branch of mathematics that extends integration and differentiation operators to noninteger 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–59} waxy crude oil,^{60,61} as well as polymers and gels,^{41,62–65} and food.^{11} 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} Considering the microscale, some observations of anomalous diffusion in complex fluids have been interpreted based on fractional viscoelastic behaviours,^{73–75} with applications to, for instance, transport movements in cell cytosol.^{76,77} Acoustic waves are also commonly used to probe or image materials. Fractional calculus appears to be a powerful tool to characterise wave propagation in complex materials exhibiting power law behaviours.^{78,79} In spite of the examples reported above, when we consider the number of powerlaw 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 powerlaw 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 springdashpot models in capturing powerlaw 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 at short and long timescales. A number of studies which originally used empirical or springdashpot 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.
2 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.^{80–83}
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.^{84} 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 compliance, 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.^{85,86} 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:

 (3) 

 (4) 
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^{86}σ(t) = G*e^{iωt} from which we can define the socalled complex modulus, or dynamic modulus:

G*(ω) = G′(ω) + iG′′(ω),  (5) 
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 modulus, are directly related. Their relationship can be expressed in the Laplace domain. For instance, by combining the Laplace transforms of eqn (3) and (4), one obtains (s)(s) = s^{−2}, where (s) and (s) are the Laplace transforms of the relaxation modulus and creep compliance respectively.^{80} Therefore, we can predict one type of behavioural mode (e.g. creep) from observing a different behavioural mode (e.g. relaxation), although such a relationship is difficult to exploit in practice due to the necessarily limited experimental time/frequency domain and the presence of noise on the data.
3 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 elementbyelement 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 Fig. 2. The creep compliance and relaxation modulus 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 compliance, 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 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.

 Fig. 2 Properties of three common viscoelastic material models.  
The differential equations relating stress and strain of traditional springdashpot models take the form of a linear ordinary differential equation with constant coefficients, and both relaxation modulus and creep compliance involve exponential functions (see Fig. 2). As an example of this, consider the stress response of a Maxwell model to a step in strain:

 (6) 
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 timescale is a useful quantity; it represents the time required for the stress to fall to 1/
e of its initial value. More importantly, the timescale 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 timescales. Some models, for example the generalised Maxwell model shown in
Fig. 3(a), give rise to an arbitrarily large number of timescales (see
ref. 80, 87 and 88 for more details about springdashpot based linear viscoelastic models).

 Fig. 3 Use of traditional viscoelastic models to analyze powerlaw responses. (a) Sketch of the Maxwell model and its generalised form with the correspondent relaxation modulus showing an exponential form. (b) Approximation of a powerlaw response σ(t) = 5·t^{−0.5} (grey line) using the generalised Maxwell model with n = 1 (blue), 2 (green), 4 (dashed purple), where for n = 1 it corresponds to the Maxwell model. Fitted model parameters in Table 2 in Appendix A. (c) Relaxation spectra of the powerlaw response and the generalised Maxwell model with n = 1, 2, and 4.  
Powerlaw viscoelastic materials correspond to systems with a broad range of relaxation/creep timescales arising from dissipation/deformation modes occurring at different time and length scales. Mathematically, powerlaw behaviour consists not of a discrete number of timescales, but a continuous distribution of timescales. For this reason, powerlaw rheology can only be approximated by the springdashpot models that yield exponential terms with discrete timescales.^{89–91} The larger the number of exponential terms, the closer the approximation.^{92–96} We demonstrate this in Fig. 3 by incrementally increasing the number of Maxwell arms in a generalised Maxwell model, each of which contributes a material timescale. With four Maxwell arms, a good approximation of a powerlaw is obtained over the time domain of the data (Fig. 3(b)). Another informative perspective on the approximation process can be gained via a relaxation spectrum, which represents the relative behavioural dominance at different times. The relaxation spectrum found here for the powerlaw relaxation modulus was approximated as^{86,97}

 (7) 
where
G(
t) is the relaxation modulus.
Fig. 3(c) shows the relaxation spectrum of the powerlaw material in logarithmic scale, which is a powerlaw distribution of timescales with the same exponent as the original relaxation data. The optimal fit of the various tested Maxwell models can be seen to yield timescales that are equally spaced in logarithmic scale. As more Maxwell arms are added, a broader coverage of the true powerlaw distribution is attained. Hence, the quality of the approximation of powerlaw responses with a generalized Maxwell model depends on both (i) the number of Maxwell components and (ii) the timewindow over which the powerlaw function is fitted.
This qualitative assessment of the above relaxation responses and their corresponding spectra shows that springdashpot models might be sufficient to approximate powerlaw 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 timescales shorter or longer than those provided in the original data.
To circumvent the above disadvantages inherent to springdashpot based approximations, empirical powerlaw expressions have been used as a simple way to capture experimental measurements.^{29,32,33,98,99} 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.
4 Fractional viscoelastic models
With regards to the modelling of powerlaw viscoelastic materials, we have now outlined the advantages and disadvantages of the phenomenological springdashpot approach and the empirical ad hoc approach. In this section we introduce an additional viscoelastic element, the springpot, which is able to capture a powerlaw 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 powerlaw signatures in various contexts.
4.1 The springpot captures powerlaw behaviour
By inspecting the powerlaw relaxation data (Fig. 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 eqn (3), the relationship between stress and strain for a powerlaw material becomes: 
 (8) 
A branch of mathematics called fractional calculus provides the definition of the generalization of the differentiation operation to noninteger order valid for a function f(t) = 0 for t < 0 as^{100} 
 (9) 
where Γ(·) is the gamma function. Redefining the coefficient A = c_{β}/Γ(1 − β) in eqn (8) yields the following result: 
 (10) 
which provides a simple constitutive equation for powerlaw 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 Fig. 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.^{101} 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.

 Fig. 4 Sketch of the fractional elementspringpot. It behaves as a spring when β = 0 and as a dashpot when β = 1.  
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.^{102,103} 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 β,^{104} and this is what we use in the current review.
4.1.1 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 powerlaw viscoelasticity and the fractional derivative has been presented before,^{105,106} here we attempt to 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. Eqn (11) can in fact be generalised to any positive real number α by recalling the extension of factorial numbers to noninteger values via the Gamma function Γ(α) = (α − 1)! 
 (12) 
If we substitute α = 1 − β into eqn (12) and , we can then rewrite the expression for the total stress reported in eqn (8) as follows 
 (13) 
Since the fundamental relation (I^{a+b}f)(t) = I^{a}(I^{b}f)(t) holds,^{100}eqn (13) can be expressed as 
 (14) 
where (I^{−β}ε)(t) is the definition of fractional derivative (D^{β}ε)(t). Equivalent to the above, we can write 
 (15) 
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 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.^{107} If we assume that the system is at rest for time t ≤ 0, the fractional derivative is given by^{108}

 (16) 
where 0 < β < 1 and C denotes a Caputo fractional derivative. 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 nonlocal 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 eqn (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.
4.1.2 Rheological behaviour of the springpot.
Despite the simple constitutive equation of the fractional springpot element, it possesses rich behavioural diversity. It is a generalization of the classical viscoelastic elements, the spring and the dashpot, and exhibits a behaviour that is intermediate between the two (see Fig. 5). By substituting a step of deformation ε(t) = ε_{0}Θ(t) (where Θ(t) is the Heaviside step function) into eqn (15) and making use of the definition of the fractional derivative in eqn (16), we can extract the relaxation modulus G(t) of the springpot:^{107} 
 (17) 
The relaxation modulus is as expected a powerlaw 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 powerlaw of exponent −β, and therefore captures very well, by design, the behaviour of powerlaw materials, with only two parameters. The creep compliance J(t) can be deduced from the expression of G(t). Both are related in the Laplace domain by (s)(s) = s^{−2}. The Laplace transform of G(t) in eqn (17) is given by (s) = c_{β}s^{β−1}. Therefore, (s) = 1/(c_{β}s^{1+β}), whose inverse Laplace gives the creep compliance^{107} 
 (18) 

 Fig. 5 Responses of a springpot when subjected to (a) constant strain and (b) a step in stress. The red curves are the responses for a spring (β = 0), the blue curves for a dashpot (β = 1) and the grey curves are for increasing values of β from 0.1 to 0.9 (orange, yellow and green curves are respectively β = 0.3, 0.5 and 0.7).  
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 β (Fig. 5(a)). When a dashpot (β = 1) is subjected to a step in strain, the stress is initially infinite but immediately dissipated afterwards (Fig. 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 (Fig. 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. Fig. 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

 (19) 
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. From eqn (19) one can notice that the recovery time scale during creep experiments depends on both the powerlaw exponent β and the time for which the load has been imposed t*. The greater the powerlaw exponent of the springpot β, the larger the time required for total strain recovery, as shown in Fig. 5(b). For example, by looking at the response of a springpot with β = 0.7 (green curve in Fig. 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 powerlaw materials. It might be useful to rewrite eqn (19) as a function of the nondimensional time τ = t/t* 
 (20) 
It can be observed that the greater the time for which the load is applied to the material, the longer is the time taken for the strain to be completely recovered. This dependency should be taken into account during the design of experiments to fully capture the rheological response of powerlaw materials.
It is also illuminating to consider the response of the springpot in the frequency domain. However, the Caputo definition for the fractional derivative is not appropriate here since the assumption that ε(t) = 0 for t < 0 does not hold anymore if we excite the springpot with a periodic deformation. To overcome this problem, we use the Generalised Liouville–Caputo formulation

 (21) 
where L denotes the generalised Liouville–Caputo fractional derivative. Note that this definition is equivalent to the Caputo fractional derivative for functions f such that f(t < 0) = 0.^{109} By now considering the oscillatory load ε(t) = e^{iωt} with eqn (21) (where ^{L}D^{β}f(t) = σ(t) and f(t) = ε(t)) the complex modulus of a springpot can be found 
 (22) 
where ω is the frequency. Separating the real and imaginary parts, we find the storage and loss moduli respectively 
 (23) 
They follow a powerlaw 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 (Fig. 6(a)). The storage and loss modulus exactly match when β = 0.5 (Fig. 6(a)). For β < 0.5 the storage modulus is always greater than the loss modulus, whilst the opposite is true for β > 0.5 (Fig. 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 , constant for all frequencies.

 Fig. 6 Storage (solid line) and loss (dash line) moduli for: (a) spring (red curve), dashpot (blue curve) and the special case of β = 0.5 (yellow curves) for which storage and loss moduli are exactly the same; (b) for β < 0.5 the storage modulus is always greater than the loss modulus (β = 0.3, orange curves), whilst the opposite is true for β > 0.5 (β = 0.7, green curves).  
4.1.3 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 ScottBlair and Veinoglou.^{44} Since then, use of the springpot has been somewhat erratic, possibly because of its nontrivial mathematical foundations. However, the springpot has found use in a diverse array of materials outside of the geological and constructionmaterials contexts; for example, in gels. An early example is the work of Winter and Chambon^{91} who derived a springpotlike 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.^{110} More recently, it has been used to capture the relaxation behaviour of human arteries^{111} and the dynamic (oscillatory) response of brain tissue for the understanding of neurological disorders.^{112} There are also several biomechanical studies which make use of a powerlaw 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 compliance.^{32,113,114} Similarly, an empirical function used to analyze the relaxation response of smooth muscle cells is equivalent to that of a single springpot.^{115}
4.2 Generalised fractional viscoelastic models
Many materials exhibit powerlaw viscoelastic behaviour. However, the powerlaw regime is often limited in time; some materials show a powerlaw response at short timescale that converges to a welldefined plateau (solid like behaviour),^{29,72} whilst others may exhibit a powerlaw response followed by a continuous flow of the material (fluid like behaviour).^{116,117} 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 springdashpot network approach can be extended to concisely capture powerlaw viscoelastic regimes. Such models are often referred to as generalised viscoelastic models and they have been mathematically characterized in the past.^{102,103,118–120} It is worth highlighting that the research field of fractional calculus is currently very active, with recent progress and ongoing work to fully characterise the viscoelastic response of fractional models in terms of novel fractional operators.^{79,121,122} 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 ref. 123 for a further review of fractional viscoelasticity in geotechnical engineering). For a comprehensive introduction to fractional calculus and fractional viscoelastic models refer to ref. 78. In the following section, we demonstrate the behavioural diversity that these generalised models are capable of.
Table 1 Examples of practical uses of generalized models
Model 
Applications 
Fractional Kelvin–Voigt 
Modelling of human prostate tissue to develop novel criteria for cancer detection.^{125} 
Modelling of the oscillatory response of canine liver.^{126} 
Modelling of the relaxation response of beast cells and tissue samples.^{40} 
Modelling of the tissue mimicking materials CF11 and gelatin.^{39} 
Modelling of the viscoelastic response of potato starch gel.^{127} 
Modelling of a threedimensional particulate gel.^{63} 
Studying of the applicability of the model to wellbore creep^{128} 

Fractional Maxwell 
Modelling of the creep behaviour of rocks.^{55} 
Modelling of the viscoelastic response of tight sandstone.^{56} 
Modelling of both colloidal^{64} and carbopol^{129} gel rheology. 
Modelling of arterial tissue.^{130} 
Modelling the properties of food gels.^{11} 
Modelling the mechanical properties of collagen gel.^{131} 

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.^{132} 
Modelling of the viscoelastic response of brain tissue.^{133} 
Modelling of the viscoelastic properties of the glassy amorphous polymer polymethylmethacrylate.^{41} 
Modelling of the mechanical properties of polymer networks with addition of fillers.^{134} 
Modelling of stress relaxation of filled and unfilled polymeric materials.^{135} 
Fractional burgers 
Viscoelastic analysis of waxy crude oil.^{60,61} 
4.2.1 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 timescale response, while the springpot with the higher exponent determines the long timescale behaviour (see Fig. 7 first column). This is true for both the relaxation and creep responses, thus for β < α < 1:

 (24) 
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 timescale (see Fig. 7 second column), thus for β < α < 1: 
 (25) 
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 liquidlike behaviour at long time scale (dominated by the dashpot, α = 1), while the Kelvin–Voigt model shows a solidlike behaviour at long timescale (dominated by the spring, β = 0).

 Fig. 7 Qualitative behaviour of two springpots in series and parallel. When two springpots are placed in series, the short time scale response is dominated by the springpot with lower exponent, while the long timescale response is controlled by the springpot with higher exponent (top figures). The contrary is true for two springpots in parallel (bottom figures). In the graphs at the bottom row, the storage modulus is the solid line, while the loss modulus is the dashed line. The parameters of the models are c_{α} = 1, α = 0.8, c_{β} = 1, β = 0.2.  
Although the above heuristics directly apply to the twospringpot models, the insights can be extrapolated to more complex models in a straightforward manner. For example, consider the standard linear solid model (Fig. 2) but with the spring k_{1} replaced by a springpot. If we look at the global network, the k_{0} spring in parallel dominates the longtime 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 timescale 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 (ESI† attached to this document), where the relaxation and creep behaviours of more complicated models are 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 that fractional models often capture the rich mechanical behaviour of powerlaw materials more accurately than springdashpot viscoelastic models while using less parameters.
4.2.2 Examples of empirical functions equivalent to fractional viscoelastic models.
In this section we show how fractional viscoelastic models can be related to a selection of widely used empirical functions. We limit the study below to functions involving timeinvariant parameters. Fractional calculus also provides a solid alternative to a number of empirical models based on nonconstant parameters, such as timedependent viscosity, to account for the complex behaviour of powerlaw materials.^{136}
4.2.2.1 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,137–142} whose complex modulus is given by: 
 (26) 
where is commonly called the hysteresivity or structural damping coefficient, β is the powerlaw 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 semifractional Kelvin–Voigt model consisting of a springpot in parallel with a dashpot. In fact, eqn (26) can be rewritten as: 
 (27) 
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}.
4.2.2.2 Two powerlaw regime.
The frequency response of single cells often exhibits two powerlaw regimes.^{38,124} A lowerpower exponent at low frequencies, followed by a higher powerlaw exponent at high frequencies. Such behaviour has been frequently been analyzed by the empirically derived superposition of two powerlaws 
G*(ω) = A(iω)^{α} + B(iω)^{β},  (28) 
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 Fig. 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 Fig. 8.

 Fig. 8 Fitting of storage modulus and loss modulus showing two powerlaw regimes using fractional Kelvin–Voigt model. Dynamic response of (a) smooth muscle cells cytoskeleton^{38} and (b) kidney epithelial cells ATPdepleted.^{124} Fitted model parameter values in Table 3.  
4.2.2.3 Powerlaw cutoff.
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 powerlaw 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 powerlaw empirical function to include an exponential regime 
G(t) = At^{−α}e^{−t/τ} + B.  (29) 
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 timescale.^{66} The novel configuration features a springpot that sustains all the load at short timescale, 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 Fig. 9(a)). The mathematical expression for the relaxation modulus of the fractional model introduced by Bonfanti et al.^{66} is 
 (30) 
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^{144,145} (see Appendix B). By fitting the fractional model to the empirical functions we can demonstrate that they are approximately equivalent (Fig. 9(a)) and therefore the empirical function can be mapped into the fractional viscoelastic framework.

 Fig. 9 Fitting of empirical functions with fractional models. (a) Biphasic empirical function describing the relaxation response of epithelial monolayers fitted with the solid fractional model for different model parameters (fitted models parameters in Table 4). (b) Modified powerlaw empirical model for the relaxation response fitted with the fractional standard linear solid model (fitted model parameters in Table 5).  
Interestingly, the same empirical model was independently developed to capture the relaxation response of one of the most widely used elastomers, polydimethylsiloxane (PDMS),^{143} and the rheological behaviour of alginatebased gels.^{146,147} The study of the liquid–solid transition during gelation has attracted significant attention in the past.^{148} Friedrich et al.^{143} reported the rheological response of stoichiometrically balanced polydimethylsiloxane (PDMS) in which the crosslinking reaction was interrupted at different times, before and after the critical point (referred to as the gel point). As shown in Fig. 10, the fourparameter fractional solid model can also be applied here to accurately captures the behaviour of the timeevolving crosslinking reaction pre and postgelation.

 Fig. 10 Storage (blue) and loss (red) moduli of crosslinked PDMS samples at different extent of reactions (t − t_{c}, where t is the current time and t_{c} the time of gelation). The xaxis is shifted by A decades to facilitate reading. The dots are the experimental data (from ref. 143), while the dashed lines represent the fitted fractional viscoelastic model shown in the inset. Fitted parameters reported in Table 7 in Appendix A.  
Physically, during gelation a polymer network forms. The critical gel point is associated with an abrupt change in viscosity,^{149} 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 A, the fluidlike behaviour of the PDMS pregel 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 (Fig. 12 in Appendix A). Postgel 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.
4.2.2.4 Modified powerlaw models.
Another empirical model, recently used for the analysis of the relaxation response of various benign and malignant cell lines, is the modified powerlaw (MPL) model^{95} which can be written as 
 (31) 
where E_{0} is the instantaneous or ‘glassy’ modulus, E_{∞} is the plateau or ‘rubbery’ modulus, τ is a time scaling, and α determines the powerlaw gradient of relaxation. In contrast to a regular powerlaw 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 (a springpot in series with a spring, all in parallel to a spring), as shown in Fig. 9(b). The relaxation modulus of the fractional standard linear solid model is defined as 
 (32) 
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.^{144,145} The similarity between the two models was also discussed by Bagley,^{150} 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 (Fig. 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 powerlaw.^{151} Interestingly, the MPL model has also been used in several asphalt and asphaltconcrete studies^{46,152,153} indicating that the fractional standard linear solid may have utility in this field.
4.2.3 Fractional models capture complex powerlaw behaviours with less parameters.
As discussed, springdashpot viscoelastic models can be used to capture powerlaw 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 powerlaw materials with only two parameters. This ability relies on the fact that while springdashpot models define specific time scales, the fractional viscoelastic element directly parameterises the relaxation spectrum, and thus the distribution of characteristic time scales. To further illustrate the ability of fractional models to accurately capture complex powerlaw responses with a minimum number of model parameters, we now reanalyse the timeresponse of different powerlaw materials, originally modelled using springdashpot 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 shorttime scale response (see Fig. 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 Fig. 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 Fig. 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 springdashpot model used in the original paper.

 Fig. 11 (a) Relaxation data from zonal articular chondrocytes^{34} fitted with a standard linear solid and fractional Kelvin–Voigt specialized by one spring. (b) Relaxation data from tomato mesocarp cells^{154} fitted with a 2 timescale generalised Maxwell model and a fractional Maxwell model specialized by one dashpot. (c) Relaxation data from PCL/bioactive glass^{155} fitted with a 2 timescale generalised Maxwell model and a fractional special model. (d) Relaxation data from collagen fibrils^{10} fitted with a 2 timescale generalised Maxwell model and a fractional special model. All fitted parameters are reported in Table 6.  
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 timescales that gives rise to the powerlaw regime at short timescale, a 2 timescale 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 (Fig. 11(b)).
The 2 timescale 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^{155} and single collagen fibrils from the extracellular matrix under tensile testing.^{10} Given that they present the same qualitative behaviour as the epithelial monolayers (powerlaw followed by exponential behaviour until a final steadystate 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 timescale responses while using one less parameter (Fig. 11(c and d)).
5 Conclusions
This review has demonstrated the ability of fractional viscoelastic models to accurately capture the rheological responses of a broad range of materials while using less parameters than traditional viscoelastic models, identifying material parameters that account for the wide distribution of timescales involved in powerlaw behaviour. Despite the limitation to the linear response, fractional models exhibit rich behaviours consistent with empirical data in both relaxation and creep experiments. 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 (ESI†) 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, RHEOS, to facilitate the fitting of experimental data and prediction of powerlaw behaviours.^{156} The software package allows nonexperts 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 powerlaw 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 powerlaw 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 powerlaw 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 powerlaw behaviours.
Conflicts of interest
There are no conflicts to declare.
Appendix A: fitted parameters
Table 2 Fitted parameters generalised Maxwell models from Fig. 3(b)
Parameters 
1 timescales 
2 timescales 
4 timescales 
k
_{1} (Pa) 
3.25 
11.50 
7.50 
η
_{1} (Pa s) 
86.57 
10.12 
1.23 
k
_{2} (Pa) 
— 
1.88 
4.71 
η
_{2} (Pa s) 
— 
103.64 
5.55 
k
_{3} (Pa) 
— 
— 
1.96 
η
_{3} (Pa s) 
— 
— 
17.14 
k
_{4} (Pa) 
— 
— 
1.06 
η
_{4} (Pa s) 
— 
— 
132.97 
Table 3 Fitted parameters fractional Kelvin–Voigt model from Fig. 8. Note that the fitted storage and loss moduli of TC7 kidney epithelial cells are normalized with respect to their values at ω = 10 Hz
Data 
c
_{
α
}

α

c
_{
β
}

β

Bovine trachea smooth muscle cells^{38} 
0.01 (Pa nm^{−1} s^{−α}) 
0.78 
0.9 (Pa nm^{−1} s^{−β}) 
0.1 
TC7 kidney epithelial cells ATPdepleted^{124} 
0.02 (s^{α}) 
0.6 
0.7 (s^{β}) 
0.07 
Table 4 Parameters of empirical function and threeelement fractional model from Fig. 9(a)

Empirical model 
Fractional model 
Curve 1 (bottom) 
A (Pa s^{β}) 
960 
c
_{
β
} (Pa s^{β}) 
1125 
α

0.15 
β

0.1 
τ (s) 
8.0 
η (Pa s) 
6650 
B (Pa) 
800 
k (Pa s) 
800 

Curve 2 
A (Pa s^{β}) 
960 
c
_{
β
} (Pa s^{β}) 
1400 
α

0.3 
β

0.21 
τ (s) 
8.0 
η (Pa s) 
6100 
B (Pa) 
1300 
k (Pa) 
1300 

Curve 3 
A (Pa s^{β}) 
960 
c
_{
β
} (Pa s^{β}) 
6000 
α

0.5 
β

0.4 
τ (s) 
8.0 
η (Pa s) 
2080 
B (Pa) 
1800 
k (Pa) 
1800 

Curve 4 (top) 
A (Pa s^{β}) 
960 
c
_{
β
} (Pa s^{β}) 
6700 
α

0.7 
β

0.54 
τ (s) 
8.0 
η (Pa s) 
3040 
B (Pa) 
2300 
k (Pa) 
2300 
Table 5 Parameters of empirical function and threeelement fractional model from Fig. 9(b)

Empirical model 
Fractional model 
Curve 1 (bottom) 
E
_{∞} (Pa) 
0.5 
c
_{
α′} (Pa s^{α}) 
0.007 
E
_{0} (Pa) 
1.0 
α′ 
0.72 
τ (s) 
1 × 10^{−3} 
k
_{1} (Pa) 
0.5 
α

0.8 
k
_{0} (Pa) 
0.5 

Curve 2 
E
_{∞} (Pa) 
0.5 
c
_{
α′} (Pa s^{α}) 
0.017 
E
_{0} (Pa) 
1.0 
α′ 
0.59 
τ (s) 
1 × 10^{−3} 
k
_{1} (Pa) 
0.5 
α

0.6 
k
_{0} (Pa) 
0.5 

Curve 3 
E
_{∞} (Pa) 
0.5 
c
_{
α′} (Pa s^{α}) 
0.051 
E
_{0} (Pa) 
1.0 
α′ 
0.42 
τ (s) 
1 × 10^{−3} 
k
_{1} (Pa) 
0.5 
α

0.4 
k
_{0} (Pa) 
0.5 

Curve 4 (top) 
E
_{∞} (Pa) 
0.5 
c
_{
α′} (Pa s^{α}) 
0.096 
E
_{0} (Pa) 
1.0 
α′ 
0.33 
τ (s) 
1 × 10^{−3} 
k
_{1} (Pa) 
0.5 
α

0.3 
k
_{0} (Pa) 
0.5 
Table 6 Fitted parameters of the traditional viscoelastic models and the fractional model from Fig. 11. Note that the fitted stress relaxation responses of PCL/bioactive glass are normalised with respect to the maximum stress
Material 
Traditional model 
Fractional model 
Zonal articular chondrocytes^{34} 
k
_{0} (nN) 
260 
c
_{
β
} (nN s^{β}) 
420 
η
_{1} (nN s) 
172 
β

0.82 
k
_{1} (nN) 
260 
k (nN) 
253 

Tomato mesocarp cells^{154} 
k
_{0} (Pa) 
8 × 10^{5} 
η (Pa s) 
6 × 10^{6} 
η
_{1} (Pa s) 
1 × 10^{6} 
c
_{
β
} (Pa s^{β}) 
1.5 × 10^{6} 
k
_{1} (Pa) 
6 × 10^{5} 
β

0.16 
η
_{2} (Pa s) 
1 × 10^{6} 


k
_{2} (Pa) 
4 × 10^{4} 



PCL/bioactive glass^{155} bottom sample 
k
_{0}

0.8 
η (s) 
1.6 
k
_{1}

0.1 
c
_{
β
} (s^{β}) 
0.15 
η
_{1} (s) 
0.01 
β

0.02 
k
_{2}

0.15 
k

0.76 
η
_{2} (s) 
1.46 



PCL/bioactive glass^{155} top sample 
k
_{0}

0.8 
η (s) 
1.5 
k
_{1}

0.1 
c
_{
β
} (s^{β}) 
0.13 
η
_{1} (s) 
0.04 
β

0.13 
k
_{2}

0.1 
k

0.8 
η
_{2} (s) 
1.02 



Collagen fibrils^{10} top sample 
k
_{0} (MPa) 
473 
η (MPa s) 
4700 
k
_{1} (MPa) 
93 
c
_{
β
} (MPa s^{β}) 
165 
η
_{1} (MPa s) 
3800 
β

0.12 
k
_{2} (MPa) 
80 
k (MPa) 
470 
η
_{2} (MPa s) 
300 



Collagen fibrils^{10} middle sample 
k
_{0} (MPa) 
414 
η (MPa s) 
5100 
k
_{1} (MPa) 
78 
c
_{
β
} (MPa s^{β}) 
162 
η
_{1} (MPa s) 
331 
β

0.21 
k
_{2} (MPa) 
62 
k (MPa) 
408 
η
_{2} (MPa s) 
2700 



Collagen fibrils^{10} bottom sample 
k
_{0} (MPa) 
395 
η (MPa s) 
4400 
k
_{1} (MPa) 
75 
c
_{
β
} (MPa s^{β}) 
120 
η
_{1} (MPa s) 
3080 
β

0.12 
k
_{2} (MPa) 
41 
k (MPa) 
390 
η
_{2} (MPa s) 
163 


Table 7 Fitted parameters of the fractional model from Fig. 10
t − t_{c} 
η (Pa) 
c
_{
β
} (Pa s^{β}) 
β

k (Pa) 
−6 
3.9 
9.4 
0.71 
0.0 
−2 
52.3 
50.5 
0.61 
0.0 
0 
336.3 
121.0 
0.54 
20.2 
2 
1122.7 
225.2 
0.46 
101.3 
6 
2963.5 
920.6 
0.36 
1923.0 

 Fig. 12 Fitted parameters of the fractional model from Fig. 10.  
Appendix B: Mittag–Leffler function
The Mittag–Leffler function arises frequently in the study of fractional viscoelasticity (and fractional calculus more generally). This appendix aims to provide an overview of the function's behavior. There are one, two and three parameter versions of the Mittag–Leffler function;^{157} only the first two are introduced here as they are most relevant. The function has a series definition and an integral definition.^{145} The series definition is more amenable to intuition and will be used here. The one parameter series representation is: 
 (33) 
The first special case of note is E_{0}(z) = 1/(1 − z), which can be confirmed by Taylor expanding the fraction about z = 0 and comparing with the series definition in eqn (33), noting that Γ(1) = 1. The second special case of note is E_{1}(z) = e^{z}, which can be seen by comparing eqn (33) to the series definition of the exponential function and noting that Γ(k + 1) = k! for k ∈ ^{0}. These two special cases illuminate the fact that the Mittag–Leffler function is capable of interpolating between power law and exponential behaviour. The two asymptotic representations are^{151} 
 (34) 
For values of 0 < α < 1 the behaviour as t → 0 approaches a stretched exponential, while for t → ∞ the Mittag–Leffler function is asymptotically equivalent to a powerlaw (Fig. 13).

 Fig. 13 Qualitative behaviour of the Mittag–Leffler function for different values of α = 0.25, 0.50, 0.75.  
The one parameter Mittag–Leffler function is in fact a special case of the two parameter function, E_{α,1} = E_{α}. The series representation of the two parameter Mittag–Leffler function is

 (35) 
As might be expected, the two parameter Mittag–Leffler function is capable of even greater behavioural diversity than the one parameter version. Although the manifold behaviours are not summarised here as they lie beyond the scope of this review, the asymptotic behaviour of the two parameter Mittag–Leffler function can be found in a recent paper.
^{158}
Acknowledgements
The authors wish to acknowledge present and past members of the Kabla and Charras labs for stimulating discussions. The authors would like to thank Prof. Manuel Duarte Ortigueira, Prof. Sverre Holm, Prof. Jordan Hristov, Dr Arran Ferandez and the referees for their comments on the preprint version of this work. AB, JLK, and AJK acknowledge the BBSRC grants BB/M002578/1, BB/K018175/1, and BB/P003184/1. JLK would like thank the George and Lillian Schiff Foundation for the PhD funding which facilitated this project. G. C. was supported by a consolidator grant from the European Research Council (MolCellTissMech, agreement 647186) and by the BBSRC grant BB/M003280.
Notes and references
 M. E. Mackay, J. Rheol., 2018, 62, 1549–1561 CrossRef CAS .
 A. Corker, H. C.H. Ng, R. J. Poole and E. GarcaTuñón, Soft Matter, 2019, 15, 1444–1456 RSC .
 R. I. Tanner, F. Qi and S.C. Dai, J. NonNewtonian Fluid Mech., 2008, 148, 33–40 CrossRef CAS .
 A. Lazaridou, D. Duta, M. Papageorgiou, N. Belc and C. G. Biliaderis, J. Food Eng., 2007, 79, 1033–1047 CrossRef CAS .
 N. Hata, A. Tobolsky and A. Bondi, J. Appl. Polym. Sci., 1968, 12, 2597–2613 CrossRef CAS .
 E. van Ruymbeke, S. Coppola, L. Balacca, S. Righi and D. Vlassopoulos, J. Rheol., 2010, 54, 507–538 CrossRef CAS .
 D. Gagnon, X. Shen and P. Arratia, EPL, 2013, 104, 14004 CrossRef .
 J. Figueroa, C. Manuel, Z. HernándezEstrada and B. RamrezWong, Cereal Chem., 2012, 89, 211–216 CrossRef CAS .
 X. Xu and J. Hou, Mech. TimeDepend. Mater., 2011, 15, 29–39 CrossRef .
 Z. L. Shen, H. Kahn, R. Ballarini and S. J. Eppell, Biophys. J., 2011, 100, 3008–3015 CrossRef CAS PubMed .
 T. Faber, A. Jaishankar and G. McKinley, Food Hydrocolloids, 2017, 62, 325–339 CrossRef CAS .
 S. Liu, L. Mo, K. Wang, Y. Xie and M. Woldekidan, Constr. Build. Mater., 2016, 105, 1–13 CrossRef CAS .
 N. Desprat, A. Guiroy and A. Asnacios, Rev. Sci. Instrum., 2006, 77, 055111 CrossRef .
 Y. Park, C. A. Best, K. Badizadegan, R. R. Dasari, M. S. Feld, T. Kuriabova, M. L. Henle, A. J. Levine and G. Popescu, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 6731–6736 CrossRef CAS PubMed .
 Y. Sun, L. G. VillaDiaz, R. H. Lam, W. Chen, P. H. Krebsbach and J. Fu, PLoS One, 2012, 7, e37178 CrossRef CAS PubMed .
 H. A. Messal, S. Alt, R. M. Ferreira, C. Gribben, V. M.Y. Wang, C. G. Cotoi, G. Salbreux and A. Behrens, Nature, 2019, 566, 126 CrossRef CAS PubMed .
 K. S. Cunningham and A. I. Gotlieb, Lab. Invest., 2005, 85, 9 CrossRef CAS PubMed .
 S. Phipps, T. Yang, F. Habib, R. Reuben and S. McNeill, Urology, 2005, 66, 447–450 CrossRef CAS PubMed .
 D. A. Moulding, E. Moeendarbary, L. Valon, J. Record, G. T. Charras and A. J. Thrasher, Blood, 2012, 120(18), 3803–3811 CrossRef CAS PubMed .
 E. S. White, Ann. Am. Thorac. Soc., 2015, 12, S30–S33 CrossRef PubMed .
 A. Fritsch, M. Höckel, T. Kiessling, K. D. Nnetu, F. Wetzel, M. Zink and J. A. Käs, Nat. Phys., 2010, 6, 730 Search PubMed .
 J. PalacioTorralba, S. Hammer, D. W. Good, S. A. McNeill, G. D. Stewart, R. L. Reuben and Y. Chen, J. Mech. Behav. Biomed. Mater., 2015, 41, 149–160 Search PubMed .
 T. W. Remmerbach, F. Wottawah, J. Dietrich, B. Lincoln, C. Wittekind and J. Guck, Cancer Res., 2009, 69, 1728–1732 CrossRef CAS PubMed .
 F. Rahimov and L. M. Kunkel, J. Cell Biol., 2013, 201, 499–510 CrossRef CAS PubMed .
 B. Chan and K. Leong, Eur. Spine J., 2008, 17, 467–479 CrossRef PubMed .
 L. Nimeskern, G. J. van Osch, R. Müller and K. S. Stok, Tissue Eng., Part B, 2013, 20, 17–27 CrossRef PubMed .
 F. Guilak, D. L. Butler, S. A. Goldstein and F. P. Baaijens, J. Biomech., 2014, 47, 1933–1940 CrossRef PubMed .
 M. M. Laronda, A. L. Rutz, S. Xiao, K. A. Whelan, F. E. Duncan, E. W. Roth, T. K. Woodruff and R. N. Shah, Nat. Commun., 2017, 8, 15261 CrossRef CAS PubMed .
 N. Khalilgharibi, J. Fouchard, N. Asadipour, R. Barrientos, M. Duda, A. Bonfanti, A. Yonis, A. Harris, P. Mosaffa and Y. Fujita,
et al.
, Nat. Phys., 2019, 15(8), 839–847 Search PubMed .
 F. Serwane, A. Mongera, P. Rowghanian, D. A. Kealhofer, A. A. Lucio, Z. M. Hockenbery and O. Campàs, Nat. Methods, 2017, 14, 181 CrossRef CAS PubMed .
 S. Nicolle, P. Vezin and J.F. Palierne, J. Biomech., 2010, 43, 927–932 CrossRef CAS PubMed .
 N. Desprat, A. Richert, J. Simeon and A. Asnacios, Biophys. J., 2005, 88, 2224–2233 CrossRef CAS PubMed .
 B. Fabry, G. N. Maksym, J. P. Butler, M. Glogauer, D. Navajas and J. J. Fredberg, Phys. Rev. Lett., 2001, 87, 148102 CrossRef CAS PubMed .
 E. Darling, S. Zauscher and F. Guilak, Osteoarthr. Cartil., 2006, 14, 571–579 CrossRef CAS PubMed .
 X. Trepat, M. Grabulosa, F. Puig, G. N. Maksym, D. Navajas and R. Farré, Am. J. Physiol.: Lung Cell. Mol. Physiol., 2004, 287, L1025–L1034 CrossRef CAS PubMed .
 A. E. Ekpenyong, G. Whyte, K. Chalut, S. Pagliara, F. Lautenschläger, C. Fiddler, S. Paschke, U. F. Keyser, E. R. Chilvers and J. Guck, PLoS One, 2012, 7, e45237 CrossRef CAS PubMed .
 E. FischerFriedrich, Y. Toyoda, C. J. Cattin, D. J. Müller, A. A. Hyman and F. Jülicher, Biophys. J., 2016, 111, 589–600 CrossRef CAS PubMed .
 L. Deng, X. Trepat, J. P. Butler, E. Millet, K. G. Morgan, D. A. Weitz and J. J. Fredberg, Nat. Mater., 2006, 5, 636 CrossRef CAS PubMed .
 F. Meral, T. Royston and R. Magin, Commun. Nonlinear Sci. Numer. Simul., 2010, 15, 939–945 CrossRef .
 H. Zhang, Q. Zhe Zhang, L. Ruan, J. Duan, M. Wan and M. F. Insana, Meas. Sci. Technol., 2018, 29, 035701 CrossRef PubMed .
 M. Alcoutlabi and J. MartinezVega, Polymer, 1998, 39, 6269–6277 CrossRef CAS .
 Y. Bouras, D. Zorica, T. M. Atanacković and Z. Vrcelj, Appl. Math. Modell., 2018, 55, 551–568 CrossRef .
 F. Barpi and S. Valente, Eng. Fract. Mech., 2002, 70, 611–623 CrossRef .
 G. S. Blair and B. Veinoglou, J. Sci. Instrum., 1944, 21, 149 CrossRef .
 G. D. Mino, G. Airey, M. Di Paola, F. P. Pinnola, G. D’Angelo and D. L. Presti, J. Civil Eng. Manag., 2016, 22, 882–889 CrossRef .
 Y. Kim, Y. Lee and H. Lee, J. Mater. Civil Eng., 1995, 7, 59–68 CrossRef CAS .
 R. Subramanian, K. Muthukumarappan and S. Gunasekaran, Int. J. Food Prop., 2006, 9, 377–393 CrossRef CAS .
 J. Lefebvre and N. Mahmoudi, J. Cereal Sci., 2007, 45, 49–58 CrossRef CAS .
 P.H. Wu, D. R.B. Aroush, A. Asnacios, W.C. Chen, M. E. Dokukin, B. L. Doss, P. DurandSmet, A. Ekpenyong, J. Guck and N. V. Guz,
et al.
, Nat. Methods, 2018, 15, 491–498 CrossRef CAS PubMed .

H. Rudolf, Applications of fractional calculus in physics, World Scientific, 2000 Search PubMed .

R. L. Magin, Fractional calculus in bioengineering, Begell House Redding, 2006 Search PubMed .
 G. S. Blair, B. Veinoglou and J. Caffyn, Proc. R. Soc. London, Ser. A, 1947, 189, 69–87 CAS .
 G. S. Blair, J. Colloid Sci., 1947, 2, 21–32 CrossRef CAS .
 H. Zhou, C. Wang, B. Han and Z. Duan, Int. J. Rock Mech. Min. Sci., 2011, 48, 116–121 CrossRef .
 F. Wu, J. F. Liu and J. Wang, Environ. Earth Sci., 2015, 73, 6965–6971 CrossRef .
 X. Ding, G. Zhang, B. Zhao and Y. Wang, Sci. Rep., 2017, 7, 11336 CrossRef PubMed .
 M. Liao, Y. Lai, E. Liu and X. Wan, Acta Geotech., 2017, 12, 377–389 CrossRef .
 C.C. Zhang, H.H. Zhu, B. Shi and B. Fatahi, Comput. Geotech., 2018, 94, 72–82 CrossRef .
 C. Ma, H.b. Zhan, W.m. Yao and H.z. Li, Water Sci. Eng., 2018, 11, 131–138 CrossRef .
 Z. Zhang, L. Hou, X. Chen, Y. Zhou, M. Liu and W. Zhou, J. Dispersion Sci. Technol., 2016, 37, 326–332 CrossRef CAS .
 L. Hou, C. Song, W. Yan and Z. Zhang, Rheol. Acta, 2014, 53, 349–356 CrossRef CAS .
 Q. Chen, B. Suki and K.N. An, J. Biomech. Eng., 2004, 126, 666–671 CrossRef PubMed .
 M. Bouzid, B. Keshavarz, M. Geri, T. Divoux, E. Del Gado and G. H. McKinley, J. Rheol., 2018, 62, 1037–1050 CrossRef CAS .
 S. Aime, L. Cipelletti and L. Ramos, J. Rheol., 2018, 62, 1429–1441 CrossRef CAS .
 J. L. Kaplan, T. A. Torode, F. B. Daher and S. A. Braybrook, bioRxiv, 2019, 565614 Search PubMed .
 A. Bonfanti, J. Fouchard, N. Khalilgharibi, G. Charras and A. Kabla, R. Soc. Open Sci., 2020, 7, 190920 CrossRef CAS PubMed .
 C. Coussot, S. Kalyanam, R. Yapp and M. F. Insana, IEEE Trans. Ultrason. Eng., 2009, 56, 715–726 Search PubMed .
 B. Carmichael, H. Babahosseini, S. Mahmoodi and M. Agah, Phys. Biol., 2015, 12, 046001 CrossRef CAS PubMed .
 Z. Dai, Y. Peng, H. A. Mansy, R. H. Sandler and T. J. Royston, Med. Eng. Phys., 2015, 37, 752–758 CrossRef PubMed .
 P. Perdikaris and G. E. Karniadakis, Ann. Biomed. Eng., 2014, 42, 1012–1023 CrossRef PubMed .
 Y. Yu, P. Perdikaris and G. E. Karniadakis, J. Comput. Phys., 2016, 323, 219–242 CrossRef PubMed .
 D. Craiem and R. L. Magin, Phys. Biol., 2010, 7, 013001 CrossRef PubMed .
 I. Goychuk, Adv. Chem. Phys., 2012, 150, 187 CrossRef CAS .
 J.H. Jeon and R. Metzler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 021147 CrossRef PubMed .
 J.H. Jeon, N. Leijnse, L. B. Oddershede and R. Metzler, New J. Phys., 2013, 15, 045011 CrossRef .
 G. Guigas, C. Kalla and M. Weiss, Biophys. J., 2007, 93, 316–323 CrossRef CAS PubMed .
 T. J. Lampo, S. Stylianidou, M. P. Backlund, P. A. Wiggins and A. J. Spakowitz, Biophys. J., 2017, 112, 532–542 CrossRef CAS PubMed .

F. Mainardi, Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models, World Scientific, 2010 Search PubMed .

S. Holm, Waves with PowerLaw Attenuation, Springer, 2019 Search PubMed .

W. N. Findley and F. A. Davis, Creep and relaxation of nonlinear viscoelastic materials, Courier Corporation, 2013 Search PubMed .
 B. Keshavarz, T. Divoux, S. Manneville and G. H. McKinley, ACS Macro Lett., 2017, 6, 663–667 CrossRef CAS .
 N. A. Bharadwaj, K. S. Schweizer and R. H. Ewoldt, J. Rheol., 2017, 61, 643–665 CrossRef CAS .
 S. Müller, M. Kästner, J. Brummund and V. Ulbricht, Comput. Mater. Sci., 2011, 50, 2938–2949 CrossRef .

T. Kailath, Linear systems, PrenticeHall Englewood Cliffs, NJ, 1980, vol. 156 Search PubMed .

D. GutierrezLemini, Engineering Viscoelasticity, Springer, 2014, pp. 23–52 Search PubMed .

R. Lakes and R. S. Lakes, Viscoelastic materials, Cambridge University Press, 2009 Search PubMed .

D. Roylance, Department of Materials Science and EngineeringMassachusetts Institute of Technology, Cambridge, MA, 2001, vol. 2139, pp. 1–37 Search PubMed .

D. R. Bland, The theory of linear viscoelasticity, Courier Dover Publications, 2016 Search PubMed .
 T. Svensson, IEEE Trans. Circuit Theory, 1973, 20, 142–144 Search PubMed .
 G. Beylkin and L. Monzón, Appl. Comput. Harmon. Anal., 2005, 19, 17–48 CrossRef .
 H. H. Winter and F. Chambon, J. Rheol., 1986, 30, 367–382 CrossRef CAS .
 N. Heymans and J.C. Bauwens, Rheol. Acta, 1994, 33, 210–219 CrossRef CAS .
 A. Halldin, M. Ander, M. Jacobsson and S. Hansson, World J. Mech., 2014, 4, 348 CrossRef .
 A. Bembey, M. Oyen, A. Bushby and A. Boyde, Philos. Mag., 2006, 86, 5691–5703 CrossRef CAS .
 Y. M. Efremov, W.H. Wang, S. D. Hardy, R. L. Geahlen and A. Raman, Sci. Rep., 2017, 7, 1541 CrossRef PubMed .

B. Wang, W. Wang, Y. Wang and L. Liu, Nano/Molecular Medicine
and Engineering (NANOMED), 2016 IEEE 10th International Conference on, 2016, pp. 51–54.

R. Christensen, Theory of viscoelasticity: an introduction, Elsevier, 2012 Search PubMed .
 Y. Kobayashi, M. Tsukune, T. Miyashita and M. G. Fujie, Phys. Rev. E, 2017, 95, 022418 CrossRef PubMed .
 N. Bonakdar, R. Gerum, M. Kuhn, M. Spörrer, A. Lippert, W. Schneider, K. E. Aifantis and B. Fabry, Nat. Mater., 2016, 15, 1090 CrossRef CAS PubMed .

K. Oldham and J. Spanier, The fractional calculus theory and applications of differentiation and integration to arbitrary order, Elsevier, 1974, vol. 111 Search PubMed .
 G. S. Blair and F. Coppen, Am. J. Psychol., 1942, 215–229 CrossRef .
 H. Schiessel, R. Metzler, A. Blumen and T. Nonnenmacher, J. Phys. A: Math. Gen., 1995, 28, 6567 CrossRef CAS .
 F. Mainardi and G. Spada, Eur. Phys. J.Spec. Top., 2011, 193, 133–160 CrossRef CAS .
 A. Jaishankar and G. H. McKinley, Proc. R. Soc. London, Ser. A, 2013, 20120284 CrossRef .

I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, Elsevier, 1998, vol. 198 Search PubMed .
 T. Surguladze, J. Math. Sci., 2002, 112, 4517–4557 CrossRef .
 M. Di Paola, A. Pirrotta and A. Valenza, Mech. Mater., 2011, 43, 799–806 CrossRef .

R. Gorenflo and F. Mainardi, Fractals and Fractional Calculus in Continuum Mechanics, 2008, vol. 1, pp. 223–276 Search PubMed .

M. Ortigueira and D. Valério, Fractional Signals and Systems, Walter de Gruyter GmbH, 2020 Search PubMed .
 B. Suki, A.L. Barabasi and K. R. Lutchen, J. Appl. Physiol., 1994, 76, 2749–2759 CrossRef CAS PubMed .
 D. Craiem, F. J. Rojo, J. M. Atienza, R. L. Armentano and G. V. Guinea, Phys. Med. Biol., 2008, 53, 4543 CrossRef PubMed .
 E. Nicolas, S. Calle, S. Nicolle, D. Mitton and J.P. Remenieras, Ultrasonics, 2018, 84, 119–125 CrossRef PubMed .
 E. Zhou, S. Quek and C. Lim, Biomech. Model. Mechanobiol., 2010, 9, 563–572 CrossRef CAS PubMed .
 P. A. Pullarkat, P. A. Fernández and A. Ott, Phys. Rep., 2007, 449, 29–53 CrossRef CAS .
 J. D. Hemmer, J. Nagatomi, S. T. Wood, A. A. Vertegel, D. Dean and M. LaBerge, J. Biomech. Eng., 2009, 131, 041001 CrossRef PubMed .
 D. Tripathi, S. Pandey and S. Das, Appl. Math. Comput., 2010, 215, 3645–3654 CrossRef .
 C. Celauro, C. Fecarotti, A. Pirrotta and A. Collop, Constr. Build. Mater., 2012, 36, 458–466 CrossRef .
 W. G. Gloeckle and T. F. Nonnenmacher, Macromolecules, 1991, 24, 6426–6434 CrossRef CAS .
 W. Glöckle and T. F. Nonnenmacher, Rheol. Acta, 1994, 33, 337–343 CrossRef .
 R. H. Pritchard and E. M. Terentjev, J. Rheol., 2017, 61, 187–203 CrossRef CAS .
 J. Y. Hristov, Front. Phys., 2018, 6, 135 CrossRef .
 J. Hristov, Eur. Phys. J. Plus, 2019, 134, 283 CrossRef .
 J. Lai, S. Mao, J. Qiu, H. Fan, Q. Zhang, Z. Hu and J. Chen, Math. Probl.
Eng., 2016, 2016, 1–15 Search PubMed .
 B. D. Hoffman, G. Massiera, K. M. Van Citters and J. C. Crocker, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 10259–10264 CrossRef CAS PubMed .
 M. Zhang, P. Nigwekar, B. Castaneda, K. Hoyt, J. V. Joseph, A. di Sant’Agnese, E. M. Messing, J. G. Strang, D. J. Rubens and K. J. Parker, Ultrasound Med. Biol., 2008, 34, 1033–1042 CrossRef PubMed .
 M. Z. Kiss, T. Varghese and T. J. Hall, Phys. Med. Biol., 2004, 49, 4207 CrossRef PubMed .
 M. D. Choudhury, S. Chandra, S. Nag, S. Das and S. Tarafdar, Colloids Surf., A, 2012, 407, 64–70 CrossRef CAS .
 P. Yu, Z. Jinzhou and L. Yongming, Petrol. Explor. Dev., 2017, 44, 1038–1044 CrossRef .
 P. Lidon, L. Villa and S. Manneville, Rheol. Acta, 2017, 56, 307–323 CrossRef CAS .
 J. J. Shen, C. G. Li, H. T. Wu and M. Kalantari, KoreaAust. Rheol. J., 2013, 25, 87–93 CrossRef .
 A. Holder, N. Badiei, K. Hawkins, C. Wright, P. Williams and D. Curtis, Soft Matter, 2018, 14, 574–580 RSC .
 M. Fraldi, A. Cugno, L. Deseri, K. Dayal and N. Pugno, J. R. Soc., Interface, 2015, 12, 20150656 CrossRef CAS PubMed .
 M. Kohandel, S. Sivaloganathan, G. Tenti and K. Darvish, Phys. Med. Biol., 2005, 50, 2799 CrossRef CAS PubMed .
 R. Metzler, W. Schick, H.G. Kilian and T. F. Nonnenmacher, J. Chem. Phys., 1995, 103, 7180–7186 CrossRef CAS .
 R. Metzler and T. F. Nonnenmacher, Int. J. Plast., 2003, 19, 941–959 CrossRef CAS .
 V. Pandey and S. Holm, Phys. Rev. E, 2016, 94, 032606 CrossRef PubMed .
 J. J. Fredberg and D. Stamenovic, J. Appl. Physiol., 1989, 67, 2408–2419 CrossRef CAS PubMed .
 J. Rother, M. BüchsenschützGöbeler, H. Nöding, S. Steltenkamp, K. Samwer and A. Janshoff, J. R. Soc., Interface, 2015, 12, 20141057 CrossRef PubMed .
 P. Cai, Y. Mizutani, M. Tsuchiya, J. M. Maloney, B. Fabry, K. J. Van Vliet and T. Okajima, Biophys. J., 2013, 105, 1093–1102 CrossRef CAS PubMed .
 B. Fabry, G. N. Maksym, J. P. Butler, M. Glogauer, D. Navajas, N. A. Taback, E. J. Millet and J. J. Fredberg, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 041914 CrossRef PubMed .
 B. A. Smith, B. Tolloczko, J. G. Martin and P. Grütter, Biophys. J., 2005, 88, 2994–3007 CrossRef CAS PubMed .
 P. Cai, R. Takahashi, K. KuribayashiShigetomi, A. Subagyo, K. Sueoka, J. M. Maloney, K. J. Van Vliet and T. Okajima, Biophys. J., 2017, 113, 671–678 CrossRef CAS PubMed .
 C. Friedrich and L. Heymann, J. Rheol., 1988, 32, 235–241 CrossRef CAS .
 H. J. Haubold, A. M. Mathai and R. K. Saxena, J. Appl. Math., 2011, 2011, 1–51 CrossRef .
 R. Gorenflo, J. Loutchko and Y. Luchko, Fract. Calc. Appl. Anal., 2002, 6(4), 491–518 Search PubMed .
 M. Nobile, V. Pirozzi, E. Somma, G. Gomez D'Ayala and P. Laurienzo, J. Polym. Sci., Part B: Polym. Phys., 2008, 46, 1167–1182 CrossRef CAS .
 M. Moresi and M. Bruno, J. Food Eng., 2007, 82, 298–309 CrossRef CAS .

H. H. Winter and M. Mours, Neutron spin echo spectroscopy viscoelasticity rheology, Springer, 1997, pp. 165–234 Search PubMed .

G. Odian, et al., Principles of polymerization, John Wiley & Sons, 2004 Search PubMed .
 R. L. Bagley, AIAA J., 1989, 27, 1412–1417 CrossRef CAS .
 F. Mainardi, Discrete Cont. Dyn. Syst., 2014, 19, 2267–2278 Search PubMed .
 H.J. Lee and Y. R. Kim, J. Eng. Mech., 1998, 124, 32–40 CrossRef .
 S. A. Forough, F. M. Nejad and A. Khodaii, Int. J. Pavement Eng., 2016, 17, 314–330 CrossRef CAS .
 Z. Li, Z. Zhang and C. Thomas, Innov. Food Sci. Emerg. Technol., 2016, 34, 44–50 CrossRef .
 A. ShahinShamsabadi, A. Hashemi, M. Tahriri, F. Bastami, M. Salehi and F. M. Abbas, Mater. Sci. Eng., C, 2018, 90, 280–288 CrossRef CAS PubMed .
 J. Kaplan, A. Bonfanti and A. Kabla, J. Open Source Softw., 2019, 4, 1700 CrossRef .
 R. Garra and R. Garrappa, Commun. Nonlinear Sci. Numer. Simul., 2018, 56, 314–329 CrossRef .
 J. Wang, Y. Zhou and D. O'Regan, Integr. Transforms Spec. Funct., 2018, 29, 81–94 CrossRef .
Footnote 
† Electronic supplementary information (ESI) available: Fractional viscoelastic models for powerlaw materials – Annex: summary of their constitutive equations, moduli & qualitative behaviours. See DOI: 10.1039/d0sm00354a 

This journal is © The Royal Society of Chemistry 2020 