Predicting muscle fatigue: a response surface approximation based on proper generalized decomposition technique

A novel technique is proposed to predict force reduction in skeletal muscle due to fatigue under the influence of electrical stimulus parameters and muscle physiological characteristics. Twelve New Zealand white rabbits were divided in four groups (n=3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=3$$\end{document}) to obtain the active force evolution of in vitro Extensor Digitorum Longus muscles for an hour of repeated contractions under different electrical stimulation patterns. Left and right muscles were tested, and a total of 24 samples were used to construct a response surface based in the proper generalized decomposition. After the response surface development, one additional rabbit was used to check the predictive potential of the technique. This multidimensional surface takes into account not only the decay of the maximum repeated peak force, but also the shape evolution of each contraction, muscle weight, electrical input signal and stimulation protocol. This new approach of the fatigue simulation challenge allows to predict, inside the multispace surface generated, the muscle response considering other stimulation patterns, different tissue weight, etc.


Introduction
During a period of activity, the force output of skeletal muscle declines. This phenomena, known as fatigue, involves multitude of processes and physiological mechanisms that remain being an object of study and analysis (Bruton et al. 2000;Allen et al. 2008). One common approach used to reproduce muscle fatigue in vitro consists of whole muscles dissected and placed in a bath solution. The in vitro protocols are very popular to study skeletal muscle fatigue without nervous system interference (Aslesen et al. 2001;Clausen and Nielsen 2007;Cairns et al. 2007Cairns et al. , 2008Goodman et al. 2009;Lunteren et al. 2011;Head et al. 2011;El-Khoury et al. 2012). Moreover, it has been demonstrated that this technique, in which muscles are stimulated electrically with all fibers being activated simultaneously and repeatedly, maintains the contractile properties of the tissue regarding calcium signaling, ion exchange, oxygen diffusion and energy metabolism (de Paula Brotto et al. 2001; Thornton et al. 2011;El-Khoury et al. 2012;Park et al. 2012;Clausen 2013a, b).
To obtain useful information in these experimental fatigue models about the rate and extent of fatigue, there have been several attempts to curve-fit the entire fatigue profile (i.e., peak maximum force/stress vs. time). For a review of different approximations in this context, the reader is referred to the work of Cairns et al. (2008). Although these fittings ensure quantitative comparisons of fatigue with different stimulation protocols, they involve averaging curves and interindividual differences are neglected (Cairns et al. 2008). Moreover, the information contained in the force development of individual contractions is not considered.
The response surface methodology (RS) has been used successfully in many research works in the biomechanics field in optimization studies (Lin et al. 2006;Nirmalanandhan et al. 2008;Eberle et al. 2013), experimental ) and computational test design (Sigal and Whyne 2010), looking for critical variables in a complex problem (Wang et al. 2005), etc. Basically, response surfaces are multidimensional surfaces fit to quantities of interest whose mathematical form allows easy interpolations to obtain another realization with a different combination of parameters. The multidimensional approximation increases its complexity and computational cost due to the wellknown "curse of dimensionality" (Ladeveze and Chamoin 2011). To reduce this order of complexity, the use of the response surface methodology combined with the technology of model reduction, known as proper generalized decomposition (PGD) (El Halabi et al. 2013b), is proposed in this work. This latter technique, based on the use of separated representations, was developed for solving multidimensional models (Ammar et al. 2006(Ammar et al. , 2007 and in the context of stochastic modeling (Nouy 2007). The technique was extended for addressing parametric models, where model parameters were considered as model extra-coordinates. This made possible the calculation of the parametric solution that can be viewed as a metamodel or a computational vademecum, to be used for real-time simulation, optimization, inverse analysis and simulation-based control (Chinesta et al. 2013).
The study reported here was conducted with the main objective of investigating the potential of the RS with the technology of PGD to predict the fatigue response of skeletal muscle. In this way, a multidimensional surface was fitted to a limited number of experimental results obtained in an in vitro animal model. Afterward, in order to validate our model, a new in vitro experiment was carried out using a different parameter combination from those used to develop the model. The results obtained were compared with the response prediction given by our model. The paper is established as follows: In the next section, the animal model and stimulation protocols to obtain the experimental data are presented. Then, a brief description of RSPGD for a general multidimensional case is presented and several numerical examples are introduced. Finally, the possibilities of the method, its cost and accuracy for different design parameters are discussed and different concluding remarks presented.

Material and methods
The experimental study was conducted on 13 male New Zealand white rabbits aged 2 months and with a body mass of 2150 ± 50 g. All experiments were approved by the University of Zaragoza Ethics Committee for the use of animals in experimentation in accordance with the provisions of the European Council (ETS 123) and the European Union (Council Directive 86/609/EEC) regarding the protection of the animals used for experimental purposes. The animals were kept in a temperature controlled room (22 ± 1 • C) with 12-h light-dark cycles and free access to water and food.

Muscle preparation
Animals were anesthetized with a medetomidine (0.14 mg/Kg), buprenorphine (0.02 mg/Kg) and ketamine (20 mg/Kg) protocol and then euthanized by an intravenous overdose of sodium pentobarbital. Immediately afterward, the animal was placed on its back and a midline incision was done to the ankle in the midsection of the knee. Tibialis Anterior (TA) was removed to access to the Extensor Digitorum Longus (EDL). The EDL muscle was carefully extracted by cutting off the distal and proximal tendons. Then, a cyanoacrylate sandpaper tab was pasted to both tendons in order to ensure a perfect attachment to the machine.

Protocol stimulation
Functional in vitro testing of rabbit muscles was carried out in a methacrylate organ bath (20 × 20 × 20 cm) designed by the authors to be installed in an electromechanical Instron Microtester 5248. The distal tendon of each muscle was strongly fixed inside the bath, and the proximal tendon of the muscles samples was also fixed to the machine actuator with a 50 N full scale load cell (Fig. 1a). Once the sample is vertically fixed, the temperature of the oxygenated Ringer's solution was maintained pumping it through a separate temperature controller and back to the organ bath. This physiological environment (27 • C and oxygen saturated solution) assured a physiological response to the electrical stimulation. To this aim, a pair of platinum plate electrodes running the length of the isolated muscle on either side was connected to a CIBERTEC CS-20 electrical signal generator. The gap between these plates can be regulated to avoid the contact with the tissue. Thus, the muscle was stimulated by the electrical field generated by the electrodes and not by direct contact (see Fig. 1a). Before the fatigue protocol, each sample was subjected to a length sweep with short active stimulation in order to determine its optimal length.
Twelve rabbits were divided in four groups (3 animals per group), and in each group left and right EDL were tested (n = 6). The amplitude of the electrical impulse was fixed to 100 V for all the groups as well as the stimulation time, fixed to 0.2 s. In the first group, EDL muscles were stimulated with a continuous signal with 100 V amplitude. The interval between train of pulses was 6 s. All the other groups received the same pulse duration (1 ms) and the same rest time (10 s). However, the frequency was 40 Hz for the second group, 60 Hz the third group and 100 Hz the fourth one. All the muscles were stimulated for 1 h. The different parameters for each group are presented in Table 1.  The amplitude of the signal was fixed at 100 V Finally, after the development of the response surface and in order to validate its predictive potential, one more rabbit was used. The signal applied to the right EDL sample consisted of 0.005 s pulses at 80 Hz during 0.2 s. The amplitude of the signal was maintained at 100 V, and the time interval between train of pulses was fixed to 10 s.
In Fig. 1b, the raw data obtained for one sample during the 1 h protocol can be observed. The acquisition software recorded force and time pair of values when the force developed by the tissue increased or decreased 1 mN.

The RSPGD approach
In this section, a brief description of the response surface using proper generalized decomposition methodology (RSPGD) (El Halabi et al. 2013b) is presented.
The RS can be easily explained as the best approximation response function f (x 1 , x 2 , . . . , x n ) partially defined by a cloud of values of that function coming from external numerical results, experiments or any other source of data. This response function is dependent on the values of n input variables or quantitative factors x = {x 1 , x 2 , . . . , x n } which are considered capable of defining a hyper surface in the bounded region Ω ⊂ R n (called the experimental region) in which we consider constrained the values of x by practical limitations.
Many types of functions have been used as approximation functions determining their associated constants using a regression analysis, like least squares technique. Other fitting methods, such as weighted least squares regression, best linear unbiased predictor used for kriging, backpropagation, mostly implemented in neural networks and many others (Nakashima 1995;Sakata et al. 2007;Arellano-Garcia et al. 2007;Park and Park 2010;Koutsourelakis 2008;Langley and Simon 1995), can be used for this proposal.
However, the fitting procedure for high-dimensional spaces is not frequently found and has not been implemented due to its intrinsic complexities (Lesh 1959). The use of high number of parameters to handle for a high-dimensional fitting increases exponentially with the number of dimensions. To avoid the so-called curse of dimensionality that appears in traditional strategies, model reduction techniques have been developed in the last years. The main characteristic of model reduction methods is that the response of complex models can be approximated by the one of a surrogate model which represents the projection of the initial model on a low-dimensional reduced basis (Nouy 2010). The difference between the different model order techniques is the way of defining and constructing this basis. The RSPGD methodology applied in this work uses a model reduction method based on separation of variables called proper generalized decomposition (PGD) (Ammar et al. 2006(Ammar et al. , 2007Chinesta et al. 2010b;Gonzalez et al. 2010). PGD has been introduced in different contexts like parameterized PDEs (Nouy et al. 2008;Nouy 2009;Nouy and Maitre 2009;Doostan and Iaccarino 2009), multiscale models (Chinesta et al. 2010a;Neron and Ladeveze 2010;El Halabi et al. 2013a), parametric modeling and structural optimization (Leygue and Verron 2010), multidimensional PDEs (Ammar et al. 2006(Ammar et al. , 2007Gonzalez et al. 2010;Nouy 2010), nonlinear models (Niroomandi et al. 2013a, b), and dynamic behavior (Gonzalez et al. 2014).
This technique is a greedy algorithm that constructs the approximation of the solution by means of a sum of products (sometimes called finite sum decomposition) of separated one-dimensional functions depending on each of the problem dimensions or parameters. Each of these functions is determined by an iterative method, with no initial assumption on their form, although usually, they are expressed as piecewise linear splines with small support as in standard linear one-dimensional finite elements.
Let ψ(x) be a scalar function of In the PGD approach, this function is approximated as: with F k i (x k ) the ith one-dimensional function of the kth variable x k that has to be computed in an implicit scheme, n the number of independent variables (dimensions), F i the product of n functions F k i and T the number of sums or terms for the approximation.
For the fitting procedure in each dimension, a least square approach is established. Given a set of pairs (x m , ψ m ) with m = 1, . . . , D known values of the function ψ for D combinations of independent parameters, find the function F(x) expressed as in (1) which minimizes E defined as: where F k i (x k ) are usually expressed in discrete form as: with N k j the standard linear one-dimensional shape function evaluated in x k , z k i j the nodal value vector of the function F k i at node j. For minimization of E, the partial derivatives with respect to the parameters of F(x) must be zero. The algorithm consists in an iterative procedure to add new terms in the finite sum until convergence of the solution. To this end, a three steps are proposed: 1. Projection step Compute the coefficients ω i from the linear system of equations derived from the minimization of the functional E with respect to the coefficients ω i The equivalent linear system may be expressed as: 2. Convergence step A convergence check for the overall solution is proposed in this stage, that is, for the already computed values of the basis functions F k i (x k ), i = 1, . . . , T k = 1, . . . , n and coefficients ω i , i = 1, . . . , T the value of the relative error must be below a certain predefined accuracy limit.
with E given by (2). If this condition is fulfilled the solving process finishes and if not, we move to the next step. 3. Enrichment step A new term T + 1 is added to the finite sum, so the new basis functions F k T +1 (x k ), k = 1, . . . , n have to be obtained. In this "enrichment stage" (Ammar et al. 2006;Gonzalez et al. 2010), the response function is then written as: Here, F k T +1 (x k ) has been substituted by R k (x k ) since the normalized basis functions F k T +1 (x k ) will be obtained by normalizing the functions R k (x k ) once the iterative process described below converges. To compute the enrichment functions, a nonlinear problem with as many equations as the number of dimensions (k = 1, 2, . . . , n) and obtained by replacing (6) into (2), must be solved.
The intervals and discretization of the variables considered in the response surface are presented in Table 2.

Experimental results
Average muscle weight, length and volume together with standard deviation are presented in Table 3. In order to obtain an estimation of the cross-sectional area (CSA) of the muscles, the technique proposed by Calvo et al. (2010) was followed, and results are shown in Table 4. All the muscles were fixed at optimal length before starting experiments in order to obtain maximal force production. Length sweep registration is shown in Fig. 2.
In Fig. 3, the experimental results are presented in the classical form in which the peak maximum stress evolution is averaged for the samples in the four groups. There were no significant differences among groups. Nevertheless, the Group 1 experimented the greatest decline, showing at t = 1000 s a reduction of around a 87.5 %. Figure 4a represents the evolution of the normalized force for the different contractions along the fatigue protocol for one of the samples (Group 1). In the figure, only every ten pulses are shown to allow a correct visualization. This evolution can be represented in a 3D plot where the third axis is the number of pulses or contractions. Plotting every ten pulses, the result is shown in Fig. 4b. Analyzing the evolution of the pulses, only the first 70 are monotonically increasing until  Fig. 2 Relationship obtained for all the muscles tested. Optimal muscle lengths were fixed to register maximum forces the impulse ends. From that point, the pulses experimented increasing sag until the 100th pulse approximately, and then, the sag remained more or less constant.

RSPGD approximation results
The experimental data were obtained using the protocol defined in Sect. 2, and a total of 12 samples were considered. The solution obtained using the RSPGD algorithm was the vectors that represent the basis functions of the approximation and its respective ω values, with which it is possible to evaluate the response at any point of the multidimensional domain.
In order to check the behavior of the response surface to obtain new values for the muscle forces, a numerical test has been developed. The response surface was evaluated at the same stimulation parameters selected for the last sample described in Sect. 2. The evolution of the maximum  Table 1. The amplitude of the signal was the same for all groups (100 V). The continuous line represents the average peak maximum value, and the gray regions are the ± standard deviation bounds . b 50th contraction obtained by the RS (blue) and experimentally (red). c Contraction evolution in a 3D representation obtained by the RS (blue) and experimentally (red) peak force for the computational and experimental results are shown in Fig. 5a. A comparison between the 50th contraction is represented in Fig. 5b. The computed relative error between the RS prediction and the experimental result is less than 13 %. Finally, the evolution of all the contractions is shown in 5c.

Discussion and conclusions
Experimental EDL force registered through the 1 h stimulation, showed a characteristic profile corresponding with fast-twitch muscles, where a significative drop of 80 % of initial force is usually observed after 1000 s of stimulation. Similar force decline was characterized before in mammals muscles with fast fiber composition (Vedsted et al. 2003;Cairns et al. 2004Cairns et al. , 2007Katz et al. 2014) which agree with fast behavior described in literature (Schiaffino and Reggiani 2011). Estimation of cross-sectional area allowed us to obtain the maximum isometric stress produced by the EDL muscles. As observed in Fig. 3, our results showed a range of tensions from around 0.4 to 0.6 MPa. In general, muscle tension in the literature reached until 0.4 MPa in mammals. Sperringer and Grange (2016) and Head et al. (2011) showed values around 0.4 MPa in mice whole muscles and the same are reported by Goodman et al. (2009) in rat EDLs. The maximum tension reached by our samples (0.6 MPa) seemed to be outside of the normal range reported by other authors. Since estimation of cross-sectional area in whole muscles is less accurate than in muscle fibers where this dimension is homogeneous through the sample, it could be possible that some of these values were overestimated.
Skeletal muscle fatigue manifests as a result of repetitive or sustained muscle contraction and can be defined as a temporary reduction in the capacity of the tissue to generate force. Many protocols have been described in the literature (Allen 2004;Brotto et al. 2004;Chan and Head 2010;Mutchler et al. 2015); however, there are several issues in both, animal and human, studies. For human, in vivo behavior study of an isolated muscle is complicated (Barbieri et al. 2014) and in vitro experiments require biopsies which are not always easy to obtain (Alghannam et al. 2014;Zampieri et al. 2015). The animal studies allow to increase the amount of samples or carry out in vitro experiments in a whole muscle (Bottinelli et al. 1991;Chan et al. 2008;Grasa et al. 2014). Nevertheless, a long-standing scientific challenge is to reduce the number of animals in experiments and replace them with reliable methods in order to obtain the same information. Computational models have been developed to predict the evolution of the force and to simulate the onset of fatigue. However, many of them only can fit the biological data obtained from the biological experiments (Tang et al. 2007;Böl et al. 2011;Grasa et al. 2014). The response surface presented here incorporates an important advantage in comparison with previous approximations because it allows to obtain new fatigue data from parameters not tested in animals. Furthermore, the mechanisms underlying muscle fatigue are not well understood yet. Due to the complexity of the phenomenon and the high number of factors included, the explanation has to take into account both the molecular actin-myosin interaction and the macroscopic phenomena observed in striated muscle as well as how conditions imposed on macroscopic scales affect actin-myosin kinetics (Minozzo et al. 2012;Röhrle et al. 2012). This response surface approximation could predict not only the force obtained in a wide range of stimulation parameters but also how muscle performance change and decline inside each stimulation, because the high precision of the acquisition system allows to register the evolution of each contraction. The results provided by this technique, when focusing on the evolution of maximum peak force, are in agreement with those obtained by other authors in muscles of different specimens (Burke et al. 1973;Darques et al. 2003;Cairns et al. 2008). Regarding the human research, this information about performance is very important in the sport field, where muscle efficiency is always the most important target (Dickerson et al. 2015;Mohr et al. 2016). Moreover, one of the main advantages of the in vitro tests developed in this work is that the central fatigue is completely eliminated (Allen et al. 2008). This could be useful when the influence of substances in the muscle performance is investigated.
The RSPGD methodology used in this work to compute the large number of parameters considered has demonstrated that few experimental data were enough to generate the approximation and to obtain reasonable good results. The use of RSPGD does not have previous constrains for the approximation, which is not the case for traditional fitting techniques where it is useful to prescribe the approximation function type to use. This characteristic gives the technique higher flexibility to fit functions with high variations. Nevertheless, the discretization must have sufficient information in each element to obtain good response, due to the local character of the technique. The principal advantage of the methodology proposed appears in high-dimensionality problems where the classical fitting procedures cannot be applied (El Halabi et al. 2013b).
The results obtained could be improved in several ways: One of them involves adding more enrichment basis functions to the approximation or adding new experimental data to enrich the sample. Another way is to change the size of the discretization for each factor. A study to optimize this discretization can be carried out, where each factor influence over the response is obtained through a previous Taguchi's design of experiments. If some factors show high variability, a better discretization for these factors must be considered, whereas if a factor shows a constant behavior, few elements will be enough.
Finally, mathematical models could lead into a better understanding of muscle fatigue to unveil the complex physiological phenomenon during continued muscle stimulation. Progress in this field, with the help of the technique proposed, would improve functional electrical stimulation (FES) techniques where muscle fatigue is currently a major drawback (Shorten et al. 2007).