Fractional Maxwell model of viscoelastic biological materials

This article focuses on fractional Maxwell model of viscoelastic materials, which are a generalization of classic Maxwell model to non-integer order derivatives. To build a fractional Maxwell model when only the noise-corrupted discrete-time measurements of the relaxation modulus are accessible for identification is a basic concern. For fitting the original measurement data an approach is suggested, which is based on approximate Scott Blair fundamental fractional non-integer models, and which means that the data are fitted by solving two dependent but simple linear least-squares problems in two separable time intervals. A complete identification algorithm is presented. The usability of the method to find the fractional Maxwell model of real biological material is shown. The parameters of the fractional Maxwell model of carrot root that approximate the experimental stress relaxation data closely are given.


Introduction
Fractional calculus is a branch of mathematical analysis that generalizes the derivative and integral of a function to non-integer order [1]. Application of fractional calculus in classical and modern physics greatly contributed to the analysis and our understanding of physico-chemical and bio-physical complex dynamical systems, since it provides excellent instruments for the description of memory and properties of various materials and processes. During the last two decades fractional calculus has been increasingly applied to mathematical modelling in physics [2,3], engineering [4,5], and especially to rheology [6,7], where fractional calculus constitute a valuable mathematical tool to handle viscoelastic aspects of systems and materials mechanics. Models involving fractional derivatives and operators have been found to better describe some real phenomena than integer-order differential equations [1][2][3][4][5], whence there are many new exciting areas of fractional models and fractional calculus applications, such as the automatic control, the modeling of biological, medical and environmental. A historical review of applications can be found in [8]. The survey on fractional models from biology and biomedicine is presented in [9], see also other papers cited therein. In recent decades fractional derivatives were found quite flexible, especially in the description of viscoelastic polymer materials [10].
Viscoelastic materials present a behaviour that implies dissipation and storage of mechanical energy. Research studies conducted during the past few decades proved that these models are also an important tool for studying the behaviour of biological materials [11]: wood, fruits, vegetables, animals tissues, etc. Viscoelasticity of the materials manifests itself in different ways, such as gradual deformation of a sample of the material under constant stress (creep behaviour), and stress relaxation in the sample when it is subjected to a constant strain. In general, viscoelasticity is a phenomenon associated with time variations in a material's response. In an attempt to describe some of the above effects mathematically several constitutive laws have been proposed which describe the stress-strain relations in terms of quantities like creep compliance, relaxation modulus, the storage and loss moduli and dynamic viscosity. Some of these constitutive laws have been developed with the aid of mechanical models consisting of combinations of springs and viscous dashpots. The classical Maxwell model, which is a viscoelastic body that stores energy like a linearized elastic spring and dissipates energy like a classical fluid dashpot is, perhaps, the most representative example of such models.
For over five decades classical exponential behaviour models have been widely applied to describe the viscoelastic properties of biological materials. Maxwell, Kelvin-Voight and Zener models are used to mathematical modelling of stress relaxation and creep processes [11][12][13]. For these models the relationship between the stress and deformation of the material is approximated though an ordinary differential or integral equations.
However, relaxation or creep processes deviating from the exponential Debye decay behaviour are often encountered in the dynamics of biological complex materials [6,13]. For such materials a stretched exponential decay KWW model (Kohlrausch-Williams-Watts) [12], hyperbolic type decay Peleg model [13] or power type behaviour models [14,15] are used to approximate experimentally obtained relaxation modulus or creep compliance data. Recently, the experimental results obtained by the authors show that the behaviour of some viscoelastic biological materials shows good agreements with that of the fractional models [9,10,16]. To this end, fractional rheological models, originally pioneered by Nutting [17] and Scott Blair [18], have proven to be a concise and elegant framework for predicting the response of complex fluids such as liquid foods using a small number of parameters [16]. By replacing the springs and dashpots of the classical viscoelastic models by the Scott Blair elements, several fractional models, including the fractional Maxwell, fractional Voigt and fractional Kelvin models, have been proposed [6,19].
In this paper fractional Maxwell model (FMM) model is considered, which relates the stress to the strain in the material by means of using differential fractional equation [3,6,7]. The FMM admit the closed form of analytical solution in terms of the known Mittag-Leffler function [20]. While fitting data to exponential sum models, like classic Maxwell model, is a very old problem, which has been studied for a long time, until now there are only a few papers concerned mainly to finding fractional Maxwell model.
We often determine the parameters in a model by obtaining the "best-possible" fit to experimental data. Common choice of the model quality measure is the mean square approximation error, leading to a leastsquares identification problem. Unfortunately, the original nonlinear least-squares task for FMM is very difficult by means of the infinite series representation of the Mittag-Leffler function [20]. Thus, in [7,21] an approach is considered, which consists on the determination of the FMM parameters by two approximate Scott Blair models identification in two separate time intervals, for small and large times respectively.
The aim of the paper is to develop a complete procedure for FMM identification using relaxation modulus data from ramp stress relaxation test.

Material
We consider a linear viscoelastic material subjected to small deformations for which the uniaxial, nonaging and izotropic stress-strain equation can be represented by a Boltzmann superposition integral [22]: where ( ) and ( ) denotes the stress and stain, respectively, and ( ) is the linear time-dependent relaxation modulus. The modulus ( ) is the stress, which is induced in the viscoelastic material described by equation (1) when the unit step strain ( ) is imposed. By assumption, the exact mathematical description of the relaxation modulus ( ) is completely unknown, but the value of ( ) can be measured with a certain accuracy for any given value of the time .

Fractional Maxwell model
Fractional Scott Blair model [7] is described by the fractional differential equation where and are the elastic modulus and relaxation time, respectively, and is non-integer positive order of fractional derivative of the strain ( ). Here, = means the derivative operator in the sense of the Caputo's fractional derivative of a function ( ) of noninteger order with respect to variable and with starting point at = 0, which is defined by [3] ( ) =  Figure 1b). To illustrate the structure of fractional models a fractional element, in addition to the standard purely elastic and purely viscous element, must be introducedsee Figure 1c. Assuming unit-step strain ( ) the uniaxial stress response of elementary fractional element (2), i.e. the time-dependent relaxation modulus ( ), is given by [6,7] Thus, the elementary fractional element (3) is uniquely described by three parameters ( , , ), as shown in Figure 1c. This approach makes it possible to include a whole range of dissipation mechanisms in a single (three parameter) rheological element. Classic viscoelastic Maxwell model is the arrangement of ideal spring in series with a dashpot (see Figure  2a) described by the first order differential equation which for unit-step strain ( ) has exponential type response with the relaxation time = ⁄ . By replacing the spring and dashpot of the classic Maxwell model by two elementary Scott Blair elements ( 1 , 1 , ) and ( 2 , 2 , )see Figure 2b we obtain fractional Maxwell model described by the fractional differential equation [6]. Without loss of generality we assume ≥ . Since, in view of (2), for the same stress ( ) the strains for these elements are described by the sum strain ( ), which is induced in the series connection of two elementary fractional elements is described by the fractional equation where the parameters of FMM (6) are functions of the parameters ( 1 , 1 , ) and ( 2 , 2 , ) of the model components given by For details of the fractional model construction see [6,7]. For the unit-step strain the solution ( ) = ( ) of FMM (6) is known for an arbitrary 1 ≥ ≥ ≥ 0 and given by the formula [6,7]: where , ( ) is the generalized two-parameter Mittag-Leffler function defined by series representation, convergent in the whole z-complex plane [3,20] , The above result is obtained in [6] by applying Fourier and Mellin transforms, while in [23] the same result is obtained by an approach involving the Laplace-Mellin transform technique. Note, that the fractional Maxwell model (7) is uniquely defined by four parameters ( , , , ), while the classic Maxwell model (4) is defined by means of only two parameters ( , ), or equivalently ( , ). A deep insight into the complex properties of FMM gives the infinite hierarchical structures of elementary fractional model (2) presented in Figure 3 composed of ideal spring and dashpot elements, which compose the classic Maxwell model.

Relaxation modulus G(t)
[Pa] derivative order are shown in Figure 4, for fixed and a few values of the respective modulus ( ) are plotted in Figure 5. Note, that the curves of FMM relaxation modulus shown in Figures 4 and 5 has the characteristic shape of the relaxation modulus of viscoelastic materials obtained in experiment. The following rule holds: the greater the parameter is, the shorter the relaxation times are for being fixed. Similar rule holds for model order , while is fixed.

Approximate fractional models
In the interpretation of the asymptotic behaviour of the FMM a fundamental role plays the asymptotic expansions of the Mittag-Leffler function. Based on (8), the relaxation modulus ( ) (7) can be rewritten as and since for small , in particular for ( ) − ≪ 1, the infinite sum of the right hand side of (9) can be truncated, the asymptotic approximation results here ≅ means 'approximately equal'. From the following asymptotic property [24] , which holds for → ∞, putting = − and = 1 − in (11) based on (7) yields whence for large times, especially for ( ⁄ ) − ≫ 1, we obtain the next asymptotic approximation of FMM in the form The two approximate rules (10), (12) are given in [6] for short times such that ≪ and for large times ≫ , respectively, based on power series representations of H-function and fact that the Mitag-Leffler function is expressible through H-function. The same asymptotic approximate models are considered in [7] and [21], but in the last paper stringent applicability conditions, derived above, are imposed. Thus, both for short as well as for large times, relaxation modulus (7) decreases almost according to the time-power elementary fractional element. The approximate models (10), (12) and the exact FMM are summarized in Figure 6, where logarithmic scale is used both for relaxation modulus and time scales. Here, the asymptotic models (10), (12), whose graph is a straight line are linear functions, which slope coefficients uniquely determine the fractional derivatives ordersfor small times and for large times, respectively.
On the basis of this observation in [7] the FMM parameters ( , , , ) of polymers PMMA and PTFE has been determined based on approximate models (10) and (12) in the time intervals defined by ≪ and ≫ , respectively. Next, in [21] the PMMA and PTFE fractional models have been corrected by applying more restrictive conditions ( ⁄ ) − ≪ 1 and ( ⁄ ) − ≫ 1. In result better fit to experimental data is obtained. In this paper a complete two-interval identification scheme based on approximate models (10) and (12) and using the linear least-squares is presented and applied to FMM of biological material identification. The scheme is inspired by [7,21] papers. However, precise applicability conditions are formulated here for the first time.

Identification
A classical manner of studying viscoelasticity is by twophase stress relaxation test, where the strain increases during the loading time interval until a predetermined strain 0 is reached at given ramp-time, after which that strain 0 is maintained constant at that value [11,12,25].
Suppose, a certain stress relaxation test performed on the specimen of the material under investigation resulted in a set of measurements of the relaxation modulus ̅ ( ) = ( ) + ( ) at the sampling instants > 0,   [12,25].
In general, identification consists of selecting within the given class of models such a model, which ensures the best fit to the measurement results. Fitting data to the original FMM (7) is a very difficult problem of nonlinear optimization, numerically difficult and often illconditioned.
Here, the linear least-squares identification routine will be used to estimate FMM parameters based on the logarithmic transformation of the experimental data and equations (10) and (12) which, respectively, yields are introduced for brevity. Denoting ( ) = ( ) and introducing the new independent variable = , equivalent linear models are obtained (the sign ≅ is neglected for simplicity) must be determined such that 1 < 2 , probably 1 ≪ 2 . The subsets 1 = {1, … , 1 } and 2 = { 2 , … , } of respective indices are chosen during the recurrent identification scheme. It is clear that the union of the sets 1 ∪ 2 does not create the set = {1, … , }. Now, classical linear least squares method can be applied to find optimal approximate models.
As a measure of the model (15) accuracy the mean sum of squares is taken the respective identification index for the second set of experimental data and log-linearized model (16) is Therefore, the least-squares identification of the loglinearized models consists of determining the model parameters minimizing the indices (18), (19) by solving the following standard optimization problems 1 (̂1 , 1 ,̂1, 1 ) = ( 1 , )∈ 2 1 ( 1 , , , 1 ), (20) 2 (̂2 , 2 ,̂2, 2 ) = ( 2 , )∈ 2 2 ( 2 , , 2 ). (21) Based on the commonly known results concerning the linear least-squares problem solution, it can be shown that the solutions to (20), (21) tasks there exist and are unique whenever the sampling instants are such that 0 < 1 < 2 < ⋯ < . The model parameters optimal in the sense of (21) are given by the known formulas: The formulas for optimal ̂1 and ̂1 , 1 are analogous. According to definitions (13) and (14) the positivity constraint can be neglected for 2 , 2 in (20), (21) optimization tasks. However, both as well as must be positive. It can be shown, based on the Czebyszev equality, that the optimal parameters ̂1 and ̂2 are positive, if and only if, the following conditions are satisfied for tasks (20) and (21), respectively.

Algorithm
Taking into account the above, the calculation of the approximate values of FMM parameters involves the following steps. 1. Perform the stress relaxation test, record and store the relaxation modulus measurements ̅ ( ), = 1, … , , corresponding to the chosen sampling instants > 0.
2. Determine the set {( , )} =1 of log-transform measurement data (17). 3. Based on the course of log-log plot of the measurement data, select in the set of measurements two separable subsets choosing 1 and 2 , 1 < 2 .
4. Compute the estimates (̂1,̂1 , 1 ) and (̂2 , 2 ,̂2) according to formulas (22), (23) and, next, compute the estimates of the relaxation time ̂1 , 2 and elastic modulus ̂1 , 2 based on (24) and (25). 5. In order to ascertain if the models with parameters (̂1,̂1 , 1 ) and (̂2 , 2 ,̂2) are a satisfactory approximation of measurement data in time intervals (0, 1 ] and [ 2 , ] compute the optimal identification indices and examine if 1 (̂1 , 1 ,̂1, 1 ) ≤ and 2 (̂2 , 2 ,̂2, 2 ) ≤ , for , a preselected small positive error. If not, change 1 and/or 2 and go to step 4. Otherwise, go to step 6. 6. In order to ascertain if the Scott Blair models (10), (12) with parameters (̂1 , 2 ,̂1 , 2 ,̂1) and (̂1 , 2 ,̂1 , 2 ,̂2) are a satisfactory approximation of the original FMM examine if If both the above conditions are satisfied, stop the procedure taking (̂1 , 2 ,̂1 , 2 ,̂2,̂1) = (̂,,̂,̂) as the FMM parameters. If not, change 1 and/or 2 and execute steps 4 and 5. The stopping rules from Step 5 guarantee the good quality of the log-linearized models in the chosen time intervals and correspond with those commonly used in the optimal identification techniques. The conditions (26), (27) are imposed to guarantee the applicability of Scott Blair models (10), (12) to approximate the original FMM. Both the pairs of conditions must be satisfied simultaneously in order to guarantee good quality of the resulted model. Note also that, in particular, from (27) the obvious necessary applicability condition of the scheme results

≫,
i.e. the experiment time must be significantly larger than the relaxation time of the FMM.
While ̂2 and ̂1 are determined independently in appropriate time intervals, the relaxation time ̂1 , 2 and the modulus ̂1 , 2 depend on the identification results in both time intervals, simultaneously.
Note, that (26) and (27) are a posteriori conditions, since the applicability of the identification procedure cannot be checked, earlier than after the experiment is performed.
The next example shows, how the identification scheme can be used for identification of FMM of real material.

Example -FMM of carrot root
A cylindrical sample of 10.9 mm diameter and 8 mm height was obtained from the root of carrot Nantejska variety [14]. During the two-phase stress relaxation test performed by Bohdziewicz [14], in the first loading phase of the time period (0,43.7) seconds the strain rate 0.000033 [ · −1 ] is used and the force induced in the specimen was recorded with the mean sampling period ∆ = 0.11007 [ ]. Next, during the second phase of the time period (43.7,1187) seconds, at constant strain the corresponding time-varying force was recorded in = 341 measurement points with sampling period ∆ = 3.34 (for details see [13,14]). Next, the relaxation modulus measurements ̅ ( ) were computed for the time interval (0,1143.3) seconds using the fast trapezoidal method of approximate relaxation modulus identification presented in [25].
Based on the log-log plot of the experiment results, the time intervals determined by 1 = 0.5 [ ] and 2 = 300 [ ] have been chosen and the estimates (̂1,̂1 , 1 ), (̂2 , 2 ,̂2) of Scott Blair models parameters, and next the estimates of the relaxation time ̂1 , 2 and elastic modulus ̂1 , 2 have been computed. Freely available [26] Matlab code for evaluating the Mittag-Leffler function has been used. The inaccurate fit of the resulting FMM to the experimental data is illustrated in Figure 7. Next, the recurrent scheme of the algorithm have been realized and successive approximations of FMM parameters have been determined.  Table 1, where both the model parameters, as well as the factors 1 , 2 (26) and Ω 1 , 2 , , defined as the mean value of the lefthand side Ω 1 , 2 of (27), are given for a few FMM approximations. Repeated composition of the scheme operations has led to parameters listed in the last row of Table 1. For these model parameters the accurate fit to experiment data is achieved, as shown in the Figure 8. The presented results convincingly prove, that the FMM give satisfactory approximation of the relaxation modulus of the material.
The last column of the table is a value of the empirical mean square model error where = ( , , , ) is the vector of FMM parameters. The model errors are not big (see Table 1), but for a few first FMM approximations the accuracy is insufficient, especially in long time region.  Next, the generalized discrete Maxwell model, which presents a relaxation of exponential type given by a finite Dirichlet-Prony series [11,12] where , and ∞ represent the elastic modulus, relaxation frequencies and equilibrium modulus (long-term modulus), respectively, has been applied to approximate the experiment data. Least-squares approximations are used, and two, three and four-parameter optimal models were determined. The parameters of the classic optimal Maxwell models and the relaxation times = 1⁄ are given in Table 2. The differences between the three Maxwell models are characterized by the ERR errors (28) summarized in Table 2, here the vector = ( , , ∞ ). Figure 9 presents experiment data and the optimal Maxwell models. It is easy to observe, that the Maxwell model is inappropriate for description of this relaxation process. Thus, the classic Maxwell models do not fully characterize the true viscoelastic behaviour of the biological material. A better fit to experimental data can be obtained, if the FMM is used. Note, that also the range of the relaxation time estimates is lower for FMM than for classic Maxwell models, even if the conditions (26)

Conclusions
In this work, a complete simple procedure is proposed for the fractional Maxwell model identification based on the discrete-time stress relaxation experimental data. The asymptotic approximate log-linearized Scott Blair models, composed of only two parameters, are used in separate time intervals. The linear least squares are applied to find the best Scott Blair models in two separable appropriately chosen time intervals. An applicability conditions of the approximate identification scheme are derived. It is shown that FMM can be used to describe the viscoelastic mechanical properties of biological materials. The effectiveness is compared with the classical two, three and four parameter Maxwell model approximations. While the idea of using the approximate Scott Blair models to find the FMM parameters is known and already appeared in [7,21], the complete identification scheme is presented in this paper for the first time.