- Email: [email protected]
Stochastic optimization for the calculation of the time dependency of the physiological demand during exercise and recovery
Computer Physics Communications 179 (2008) 888–894
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Stochastic optimization for the calculation of the time dependency of the physiological demand during exercise and recovery M.S. Zakynthinaki a,b,∗ , J.R. Stirling b a b
Instituto de Ciencias Matemáticas, CSIC - UAM - UC3M - UCM, c/Serrano 121, 28006 Madrid, Spain Facultad de Ciencias de la Actividad Física y del Deporte, Universidad Politécnica de Madrid, Avd. Martin Fierro s/n, 28040 Madrid, Spain
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 25 January 2008 Received in revised form 3 June 2008 Accepted 31 July 2008 Available online 5 August 2008 PACS: 05.45.Tp 07.05.Tp 47.10.Fg 47.63.Cb 87.19.Hh
The stochastic optimization method ALOPEX IV is successfully applied to the problem of estimating the time dependency of the physiological demand in response to exercise. This is a fundamental and unsolved problem in the area of exercise physiology, where the lack of appropriate tools and techniques forces the assumption and the use of a constant demand during exercise. By the use of an appropriate partition of the physiological time series and by means of stochastic optimization, the time dependency of the physiological demand during heavy intensity exercise and its subsequent recovery is, for the first time, revealed. © 2008 Elsevier B.V. All rights reserved.
Keywords: ALOPEX stochastic optimization Demand Modelling Physiological time series Heart rate Nonlinear dynamical systems
1. Introduction 1.1. A brief introduction to ALOPEX IV The version IV of the stochastic optimization method ALOPEX (ALgorithm Of Pattern EXtraction) has been proven to be the fastest and easiest in its implementation [1–3] when compared to other, previous versions of the same method (versions ALOPEX I, ALOPEX II [3–5] and ALOPEX III [3,6,7]). All ALOPEX methods, however, are shown to be very powerful, especially in optimization problems of many variables in real time. They are fast, effective and very easy to implement. The main advantage of ALOPEX methods is that no knowledge of the dynamics of the system or of the functional dependence of the cost function on the control variables, is required. To give a brief introduction to ALOPEX IV stochastic optimization, let us assume that the aim is to find the values of the N control variables x1 , x2 , . . . , x N that maximize function
f (x1 , x2 , . . . , x N ) = f (x) (the cost function). The method of ALOPEX IV works as follows: if xi (k) is the value of the ith control variable after the kth iteration and f (k) (x1 (k) , x2 (k) , . . . , x N (k) ) is the value of the control function after the kth iteration, then, the value of the ith control variable in the next (k + 1)th iteration is (k+1)
xi
(k)
= xi + c xi (k)
f (k) (k) + gi , | f (k−1) |
i = 1, . . . , N
(1)
where
xi ≡ xi (k) − xi (k−1) , (k)
f (k) ≡ f (k) − f (k−1) , and
f (k−1) ≡ f (k−1) − f (k−2) .
*
Corresponding author at: Instituto de Ciencias Matemáticas, CSIC - UAM - UC3M - UCM, c/Serrano 121, 28006 Madrid, Spain. Tel.: +34 626575554. E-mail addresses: [email protected] (M.S. Zakynthinaki), [email protected] (J.R. Stirling). 0010-4655/$ – see front matter doi:10.1016/j.cpc.2008.07.012
©
2008 Elsevier B.V. All rights reserved.
A proof of convergence of ALOPEX IV, in the absence of noise, can be found in [8]. For more detailed discussion on ALOPEX IV see [3,8].
M.S. Zakynthinaki, J.R. Stirling / Computer Physics Communications 179 (2008) 888–894
1.2. Physiological background Let us denote a physiological variable as s( v , t ) (referring to either the heart rate, or the rate of change of Oxygen uptake), assuming that the current value of s depends on the intensity of the particular exercise (velocity v) and the time passed after the beginning of the exercise, t. The understanding of s( v , t ) in response to exercise, is a fundamental area of exercise physiology and has many applications in sport, see [9–13] as well as medicine and health in general [14]. The values of s( v , t ) always lie within the physiological limits of that particular variable, i.e. smin s( v , t ) smax (for example there is a maximum rate that the heart is able to sustain, as well as a minimum rate, which corresponds to absolute resting values). In the present study we focus on exercises of constant intensity (v constant), as this simplest case is currently the main area of research in exercise physiology. The term ‘on-transient’ kinetics refers to the change in the values of s( v , t ) in response to an increase in the intensity of the exercise, while the term ‘off-transient’ refers to the kinetics of s( v , t ) as the body recovers to a new lower value of s, following a decrease in the exercise intensity. It has been observed that the on transient kinetics of s( v , t ) are slightly different for different constant exercise intensities. For moderate intensities, s( v , t ) follows an exponential-like rise [15,16] until it reaches a steady state. For very heavy intensity exercise the final steady state value is the maximum value smax ; there is, however, a delay (a slowing down of the kinetics) in reaching this value. For severe intensity exercise the values of s( v , t ) rise very steeply and approximately exponentially, until they are limited by the maximum value smax , which then becomes the steady state for the remaining time the exercise can be carried out at. Note that if the exercise is too severe then the subject may have to stop due to fatigue before they reach the smax . Let us assume the velocity and time-dependent function D ( v , t ) 0 that describes the physiological demand for that particular exercise (see also [8,17,18]). It is generally assumed that, for constant intensity exercises where the demand is such that D ( v , t ) smax , the value of the demand is constant and equal to the steady state value that s( v , t ) finally reaches [19] (the asymptote of the physiological time series). This time-independent value D ( v ) depends on the intensity of the particular exercise can easily be obtained from the time series of s( v , t ). For severe or very high intensity exercise there is D ( v , t ) smax . In the literature it is commonly assumed that, even in these cases, the demand is only a function of v and does not depend on time [26]. However this assumption can be shown to be an approximation which is not valid in general [8]. Indeed, for higher values of demand, where D ( v , t ) > smax , the function of demand is probably also a function of time. It is known that the energy cost of accelerating from a particular value of velocity to another value of velocity is greater, due to inertia, than the energy cost of remaining at that velocity. It is also known that the efficiency of muscle decreases with time, during constant intensity exercise [20–24]. For prolonged exercise and at a given constant intensity the energy cost increases due to a loss of efficiency during the exercise. Efficiency has also been shown to be dependent on exercise intensity, with the mechanical efficiency being lower in exhaustive exercise than in sub-maximal exercise [19,25,26]. The on-transient demand should be therefore considered to be a function of time, both in the initial and also the final stages of a bout of heavy intensity exercise. In the case of recovery (off-transient), when the intensity of the exercise is reduced to a new constant level, there is also a time dependency in the demand. It has been observed that in the initial stages of recovery following a very heavy work load, high values of s( v , t ) [27] persist for some seconds before they rapidly drop (see
889
also [28]). The reduction in the heart rate then begins to decrease as we approach the new recovery demand. There is then an ultraslow reduction in the heart rate depending on the severity and duration of the previous exercise; this can last for several hours [9]. It is known that the time to reach resting level of heart rate or oxygen uptake can be very long following very intense exercise, in fact it may not reach resting levels the same day of the exercise [9, 26]. These phenomena suggest that the demand remains elevated for a period of time during recovery. 1.3. The present study A previous study by the authors [8] showed that the ontransient heart rate demand (and also the oxygen demand [29]) can be modelled, using the model of Stirling et al. [17,18], as a constant, only for sufficiently low exercise intensities and for the off-transients following sufficiently low on-transient exercise intensities. The same study showed, however, that these assumptions break down for higher exercise intensities; a model curve assuming constant, time-independent demand could not be optimally fit to the time series of s( v , t ) for heavy intensity exercise [8]. This result justified the expectation for the demand stated in the paragraphs above. The demand, therefore, during exercise and recovery for physiological variables such as the heart rate and oxygen uptake is, in general, a function of time and not a constant, as is usually assumed. The present study applies ALOPEX IV stochastic optimization to the problem of optimizing the fit of the model [17, 18] to physiological time series data (see also [8]) with the aim of revealing, through the process of optimization, the time dependent nature of the physiological demand during exercise and recovery. For the purposes of our study heart rate time series data corresponding to heavy intensity exercise and its subsequent recovery is used. The data consists of two sets of exercise and recovery for two different heavy exercise intensities. This data is the same data that was used in [8] to show that the assumption of a constant demand during heavy intensity exercise and the recovery that follows is incorrect. Our method assumes a partition of the physiological time series in respect to time and calculates a value of the demand for each subset of the data (see Section 3). 2. The optimization problem 2.1. Modelling the kinetics of s( v , t ) The kinetics of a physiological variable s( v , t ) have been recently modelled [8,17,18] using a set of coupled ordinary differential equations of the form
˜ s( v , t ) − smin s˙ ( v , t ) = A v˙ = I (t )
B
C
smax − s( v , t )
E
D ( v , t ) − s( v , t )
(2) (3)
˜ , B , C and E where v is the velocity of the particular exercise, A are parameters that characterize the subject’s current fitness condition (for a discussion on the meaning of these parameters see [8, 17,18]), D ( v , t ) is the demand of the particular exercise (assumed to be dependent on the velocity of the exercise and time) and smin s( v , t ) smax . The function I (t ) describes the rate of change of velocity; in the present study we will assume that I (t ) = 0 and therefore v = constant for a particular exercise. For numerical reasons and without loss of generality, we assume the dimensionless normalized variable s˜ ( v , t ) ≡
s( v , t ) − smin smax − smin
890
M.S. Zakynthinaki, J.R. Stirling / Computer Physics Communications 179 (2008) 888–894
so that 0 s˜ ( v , t ) 1, and a function of a normalized demand
˜ (v , t) ≡ D
that, in accordance with the analysis of the model presented in [18], during the modelling and optimization process of the present work:
D ( v , t ) − smin smax − smin
see also [8,29]. The dynamics of the dimensionless normalized variable s˜ ( v , t ) will then be described by the equation
B
s˙˜ ( v , t ) = A s˜ ( v , t )
C
1 − s˜ ( v , t )
E
˜ ( v , t ) − s˜ ( v , t ) D
(4)
˜ (smax − smin ) B +C + E −1 . where A is the dimensionless parameter A = A Assuming that the parameters A , B , C and E of Eq. (4) are ˜ ( v , t ) that describes the demand, known and so is the function D Eq. (4) can be integrated in respect to time to give the time series {˜sm } N −1 that model the response of the physiological variable k k=0 s˜ (normalized heart rate), see also [8,18,29]. 2.2. Fitting the model to physiological time series data—The optimization parameters
• the value of the parameter E is equal to E = 1 • the value of the parameter C lies within the limits cos(C π ) > 0. Taking the above into account, it is easy to see that the shape of }N −1 that models the data depends on the values of the curve {sm k k=0 the parameters A , B and C of Eq. (4). As mentioned before these values are characteristic of the particular subject, do not change with exercise intensity and give an indication of the subject’s level of fitness, see also [8,17,18,29]. When the values of A , B and C are kept constant for a particular subject, what is left to control the shape of the model curve are the values of D ( v , t ) which, in the present study, are assumed to be time-dependent. For an exercise of a particular constant intensity the control variables of the optimization are the values of the function D (t ).
In the following sections we drop the tildes for simplification. This way the variable s( v , t ) will refer, from now on, to the normalized variable s˜ ( v , t ) and the function D ( v , t ) will refer to the ˜ ( v , t ). normalized demand D
3. The method
Let us assume that the time series describe the basic response pattern of the data (the smooth curve describing the kinetics of the data [30]). For each point k = 0, . . . , N − 1 the sum of the vertical distances (residuals R) between the smooth curve describing the data and the curve {sm } N −1 provided by the model is k k=0
As the aim of the present work is to emphasize the need of time dependent demand, we will use the same data sets used in our previous study [8] that correspond to heart rate data. We note that the experiment was carried out on a tartan track and followed to as good an approximation as possible a square wave protocol consisting of five work periods, of four laps each (for more details on the experimental protocol see [8]). The velocities of the on on sets were, v on 1 = 13.4 Km/h, v 2 = 14.4 Km/h, v 3 = 15.7 Km/h, on v on = 17 . 0 Km / h and v = 17 . 9 Km / h, where the superscript on 4 5 refers to the on-transient exercises. Every exercise was followed by a static recovery of 10 minutes, during which the subject lay hor-
f N −1 {sk }k=0
R=
N −1
f
sk − sm k
2
.
k=0
As described in [8,17], the curve that the model calculates is considered to be an optimal fit to the smoothed data curve {sk }kN=−01 if the value of R is minimum. This is an optimization problem equivalent to the problem [8] of maximizing a cost function f defined as f
f ≡ −R = −
N −1
f
sk − sm k
2
.
(5)
k=0
The aim of the optimally fitting the model curve to the data is achieved via correctly estimating the demand that corresponds to the particular data set and appropriately modifying the shape of the curve provided by the model, i.e. by changing the values of the parameters A , B , C and E of the model [17,18]. We note
3.1. The data
off
izontally and still on the floor (v i = 0, for i = 1, . . . , 5 where the subscript off refers to the off-transient recovery period). In what follows we focus on the study of the heart rate data of exercises 4 and 5, as this data corresponds to heavy intensity exercise. Indicatively, Fig. 1 shows the time series data (raw and normalized) of exercise 4 for the total time of the exercise (338 s on-transient and 600 s off-transient). In this figure, the heart rate data of 20 s before the start of the on-transient part is also shown. Fig. 2 presents, as an indication of the speed time series data, the actual time series of the velocity during the on-transient period of exercise 4. As can be seen, the velocity is as close to
Fig. 1. Raw normalized data, exercise 4. On-transient: 321 s, off-transient: 600 s.
M.S. Zakynthinaki, J.R. Stirling / Computer Physics Communications 179 (2008) 888–894
891
Fig. 2. Exercise 4, on-transient. Velocity time series.
constant as possible in an experimental situation, including measurement error.
the length of the subsets will be too large to consider a constant demand • an M too large (corresponding to a large number of small subsets) will result in miscalculation of D i (tk ): the length of the j subsets will not be enough to correctly calculate a value of D i .
3.2. The cost function The model of Eq. (4) was optimally fit to the data sets described in Section 3.1 above, as described in our previous study [8]. The assumption of a constant, time-independent demand resulted in an optimal fit of the model only for the data sets 1–3, the ones that corresponded to medium exercise intensity. The optimal parameter values, for the particular subject were found to lie within a narrow parameter neighborhood. More specifically, the best parameter combination that optimally fit the model to the data was found to be A = 0.54,
B = 1.63,
C = 1.75.
With these parameter values Eq. (4) becomes
1.63
s˙ ( v , t ) = 0.54 s( v , t )
1.75
1 − s( v , t )
D ( v , t ) − s( v , t ) .
(6)
Let us consider the ith filtered data set [30], where i = 4, 5 as we consider only the data of exercises 4 and 5. This data consists of the N points {si k }kN=−01 and the time takes the discrete values tk , k = 0, . . . , N − 1. The differential, in respect to time, Eq. (6) has the demand D i (tk ) as a variable: f
f 1.63
si˙m = 0.54 si k k
f 1.75
1 − si k
f
D i (tk ) − si k
k = 0, . . . , N − 1, i = 4, 5.
(7)
{s i m } N −1 k k=0
Time integration of (7) yields the time series that provide the model curve. We note here that for the numerical integration, the fifth-order Cash–Karp Runge–Kutta method [31] was used. The time series of the demand D i (tk ) are calculated by the maximization of the cost function as described by Eq. (5). 3.3. The partition of the data sets Let us consider a partition of the data set into M parts. The partition is done in respect to time, such that the time duration (t 1 −t 0 ) . For the jth subset of the data, where of each subset is N −M
By trial and error it was found that the optimal partition was such that the time duration of each data subset was within the interval [18 s, 40 s]. 4. Results 4.1. Estimating the values of the demand from the time series data In the sections that follow we indicatively present the results of a time partition of 20 s. We concentrate on three main areas of the heart rate data: the first part of the on-transient, the very last part of the on-transient and the initial stages of the off-transient and finally, the rest of off-transient time series. As mentioned in Section 3.1, Fig. 1 indicatively shows the time series data of exercise 4, both for the on- and the off-transient. In Fig. 1 we have divided the data into 5 sections: (1) section A includes the first 200 s (approximately) of the ontransient (2) section B includes the second part of the on-transient, from approximately 200 s up to 20 s before the end of the ontransient (3) section C includes the last 20 s of the on-transient (4) section D includes the first 20 s of the off-transient (5) section E includes the remaining of the off-transient data. In Fig. 1 the heart rate data for t ∈ [−20, 0] before the start of the exercise are also shown. It should be noted here that the on-transient part of both exercises 4 and 5 was preceded by a rest period for which it was possible to estimate a constant value of the demand. This value was equal to the asymptotic value of s during rest, which was found to be equal to 0.227 for both exercises 4 and 5.
j
j = 1, . . . , M, we consider a demand of constant value D i that will be calculated via the optimization process described in Section 3.2. The time series D i (tk ) of the demand of the ith data set will then j be composed of all the D i values, j = 1, . . . , M. It should be noted here that the process of finding the demand D i (tk ) is strongly dependent on the choice of M, as
• an M too small (corresponding to a small number of large subsets) will result in loss of detail in the time series of D i (tk ):
4.1.1. Estimating the demand during the first part of the on-transient (area A, Fig. 1) Eq. (8) gives the estimated values of the demand, for the time intervals of the on-transient part of exercise 4,
j on D4
=
⎧ 1.0723, ⎪ ⎪ ⎨ 0.8951, ⎪ 0.9132, ⎪ ⎩ 0.9481,
j=1 j=2 j=3 j4
time partition = 20 s
(8)
892
M.S. Zakynthinaki, J.R. Stirling / Computer Physics Communications 179 (2008) 888–894
and Eq. (9) shows the results for the on-transient part of exercise 5:
j on
D5
=
⎧ 1.1102, ⎪ ⎪ ⎨ 0.9362, ⎪ 0.9435, ⎪ ⎩ 0.9591,
j=1 j=2 j = 3, 4
time partition = 20 s.
(9)
j5
The results were as expected: as can be seen by observation of Eqs. (8) and (9), the demand reaches an highly elevated value straight after the beginning of the exercise. This value drops after about 20 s after the beginning of the exercise and then rises again. 4.1.2. Estimating the demand during recovery (area E, Fig. 1) The estimated values of the demand calculated for the offtransient part of exercise 4 and for a time partition of 20 s are given in Eq. (10) below. It should be noted here that, in Eq. (10), there is j = 1, . . . , 30, since the total time of the off-transient part of each exercise was 600 s. It should also be noted that the values of the demand for the first 20 s of the off-transient are found by a combination of the results of the present section and the results of the following Section 4.1.3.
⎧ 0.2901, j = 1, 2 ⎪ ⎪ ⎪ ⎪ . 2852 , j=3 0 ⎪ ⎪ ⎪ ⎪ ⎪ j = 4, 5, 6 0.2651, j off ⎨ D4 = 0.2453, j = 7, 8, 9, 10 ⎪ ⎪ ⎪ . 2335 , j = 11, 12, 13, 14, 15 0 ⎪ ⎪ ⎪ ⎪ ⎪ 0.2208, j = 16, 17, 18, 19 ⎪ ⎩ j 20 0.2134, time partition = 20 s.
(10)
Eq. (11) gives the estimated values of the demand as calculated for the off-transient part of exercise 5:
j off D5
=
⎧ 0.4294, ⎪ ⎪ ⎪ ⎪ 0.3652, ⎪ ⎪ ⎪ ⎪ 0.3064, ⎪ ⎪ ⎪ ⎨ 0.2827,
j=1 j=2 j=3
j = 4, 5 ⎪ . 2645 , j = 6, 7 0 ⎪ ⎪ ⎪ ⎪ . 2538 , j=8 0 ⎪ ⎪ ⎪ ⎪ ⎪ ≈ 0.2500, 9 j 22 ⎪ ⎩ ≈ 0.2450, j 23
time partition = 20 s.
(11)
As was expected (see Section 1.2) the heart rate demand is found to be elevated at the beginning of the off-transient period and slowly decreasing until it reaches an almost steady state (if the recovery period could last for longer, the values of the demand would asymptotically approach the minimum rest value of D = 0). 4.1.3. Estimating the demand during the last part of the on-transient and the initial stages of the off-transient (areas C and D, Fig. 1) In cases where the exercise intensity is heavy (as is the case for exercises 4 and 5) an estimation of the values of the demand for the last part of the exercise is not possible. It could be said that the last part of the on-transient time series data, where the values of s( v , t ) reach an asymptote (which is near smax , as the intensity is heavy) is a “blind part” of the time series data. In the data of exercise 4 (see Fig. 1) this “blind part” could be said to start after about 200 s from the start of the exercise. Our method estimates the value of the demand at the very last part of the on-transient time series data by making use of the fact that the kinetics of s( v , t ) can be described as a dynamical system, see Eq. (4). When time is allowed to go backwards this system will be therefore governed by the same dynamics. The kinetics backwards in time of the first part of the off-transient time
series will show what demand preceded the off-transient and will this way reveal the heart rate demand at the very last part of the on-transient exercise. We followed the same procedure described in Sections 4.1.1 and 4.1.2; the only difference is that we allowed time to go backwards, i.e. the procedure started from the off-transient and finished at the on-transient. Using the off-transient data of the first 20 s of the off-transient part (area D, Fig. 1) a value for the demand was estimated for the last 20 s of the on-transient time series data (area C , Fig. 1). The estimated (normalized) values of the demand for the last part of the on-transient of exercises 4 and 5 were found to be equal to 1.04 and 1.09 respectively. Note that these values are also the values of the demand for the beginning of the off-transient: during the first 20 s of the off-transient, the demand drops from 1.04 to 0.2901 and from 1.09 to 0.4294 for exercises 4 and 5 respectively (see Eqs. (10) and (11)). 4.2. The time dependency of the demand Figs. 3 and 4 summarize the results presented in the previous Sections 4.1.2, 4.1.1 and 4.1.3. The dashed vertical lines in these figures denote the end of the on-transient and the beginning of the recovery (off-transient) period. A smooth curve has been fitted through the estimated values of the demand, both for the on- and the off-transient part of the exercises. As discussed before, see Section 4.1.3, no information is provided regarding the “blind part” of the time series data (area B, Fig. 1). It is assumed, however, that the values of the demand throughout the exercise follow a smooth increase, as the exercise intensity was constant. In Figs. 3 and 4 we have estimated the time dependency of the demand for this “blind” area by smoothly connecting the estimated values of the demand for the first t 200 s part of the on-transient (area A) with the estimated elevated value of the very last part of the on-transient (area C ), see the dashed curves in Figs. 3 and 4. Indicatively we present in Fig. 5 the (excellent) fit of the model curve on the on-transient time series data of exercise 4. The values of the demand used for the model curve were the ones given in Eq. (8). We also present in Fig. 6 the fit of the model to the off-transient time series of exercise 5 using the values of the demand of Eq. (11). For comparison, Fig. 7 presents the fit of the model when a demand of constant value ( D 5 )off = 0.25 is considered throughout the off-transient period. 5. Conclusions We present here a novel method for estimating the time dependency of the demand during exercise and recovery, which is based on stochastic optimization. In particular, the stochastic optimization method ALOPEX IV was successfully applied to the problem of optimally estimating the values of the demand of a set of heart rate data. The cost function of the optimization was a measure of the optimal fit of a recently developed dynamical systems model [18] on the time series of the data. For this reason, the data was partitioned in respect to time and the value of the demand, for each time interval, was treated as an optimization variable. The same data set was used in a previous study by the authors [8] for the calculation of the optimal values of the model’s parameters. For the study presented in [8] a constant value for the demand was assumed during exercises of constant intensity. The results, however, showed that such an assumption is not valid for exercises of heavy intensity. The study therefore concluded that a time dependent demand should be considered for optimally modelling the kinetics of the physiological time series data (in this case, heart rate data). As can be observed in Figs. 5 and 6 (and
M.S. Zakynthinaki, J.R. Stirling / Computer Physics Communications 179 (2008) 888–894
Fig. 3. Exercise 4. Estimated normalized demand. On-transient: 338 s, off-transient: 600 s.
Fig. 4. Exercise 5. Estimated normalized demand. On-transient: 321 s, off-transient: 600 s.
Fig. 5. Exercise 4. 338 s on-transient heart rate data (normalized) and model curve. Time dependent demand.
893
894
M.S. Zakynthinaki, J.R. Stirling / Computer Physics Communications 179 (2008) 888–894
Fig. 6. Exercise 5. 600 s off-transient heart rate data (normalized) and model curve. Time dependent demand.
Fig. 7. Exercise 5. 600 s off-transient heart rate data (normalized) and model curve. Constant demand.
especially in Fig. 7 where a constant value of the demand is used to model the off-transient kinetics of the exercise), the success of the method we present here is obvious. The novelty of our method is summarized in Figs. 3 and 4 that show the time dependency of the demand for the data of exercises 4 and 5 respectively. The current literature in exercise physiology assumes a constant demand for an exercise of constant intensity, even though there are many errors in this assumption. The present work shows however that, by the use of a simple method such the one we propose here, it is possible to reveal the values of the demand that are hidden in the kinetics of the data. This way, no assumptions on the values of the demand are necessary and therefore new physiological observations or solutions to classical problems in exercise physiology are now possible. Acknowledgements This work was supported by the programs Ramón y Cajal 2004 and I3 2006, Ministerio de Educación y Ciencia, Spain. References [1] M.S. Zakynthinaki, Y.G. Saridakis, Numer. Algorithms 33 (1) (2003) 509. [2] M.S. Zakynthinaki, Y.G. Saridakis, Comp. Phys. Commun. 150 (3) (2003) 274. [3] M.S. Zakynthinaki, Stochastic optimization for adaptive correction of atmospheric distortion in astronomical observation, PhD thesis, Technical University of Crete, 2001. [4] E. Harth, E. Tzanakou, Vision Research 14 (1974) 1475. [5] T. Tzanakou, R. Michalak, E. Harth, Biol. Cybern. 35 (1979) 161. [6] T. Kalogeropoulos, Y.G. Saridakis, M.S. Zakynthinaki, Comp. Phys. Commun. 99 (1997) 255. [7] Y.G. Saridakis, M.S. Zakynthinaki, T. Kalogeropoulos, Internat. J. Appl. Sci. Comp. 5 (3) (1999) 252. [8] M.S. Zakynthinaki, J.R. Stirling, Comp. Phys. Comm. 176 (2) (2007) 98. [9] P.O. Astrand, K. Rodahl, H.A. Dahl, S.B. Stromme, Textbook of Work Physiology: Physiological of Bases of Exercise, Human Kinetics, 2003.
[10] G.A. Brooks, T.D. Fahey, T.P. White, K.M. Baldwin, Exercise Physiology, Human Bioenergetics and its Applications, third ed., McGraw-Hill, 1999. [11] D.E. Martin, P.N. Coe, Better Training for Distance Runners, Human Kinetics, Champaign IL, 1997. [12] T. Noakes, Lore of Running. Discover the Science and Spirit of Running, fourth ed., Human Kinetics, Champaign IL, 2002. [13] D.C. Poole, A. Jones, Oxygen Uptake Kinetics in Sport, Exercise and Medicine: Research and Practical Applications, Routledge, 2004. [14] K. Wasserman, J.E. Hansen, D.Y. Sue, R. Casaburi, B.J. Whipp, Principles of Exercise Testing and Interpretation, Including Pathophysiology and Clinical Applications, third ed., Lippincott Williams and Wilkins, 1999. [15] A.V. Hill, L. Lupton, Q. J. Med. 16 (1923) 135. [16] R. Margaria, F. Mangili, F. Cuttica, P. Cerretelli, Ergonomics 8 (1965) 49. [17] J.R. Stirling, M.S. Zakynthinaki, I. Refoyo, J. Sampedro, J. Nonlin. Math. Phys. (2008), in press. [18] J.R. Stirling, M.S. Zakynthinaki, B. Saltin, Bul. Math. Biol. 67 (5) (2005) 989. [19] B. Saltin, in: S. Maehlum, S. Nilsson, P. Renstrom, (Eds.), An Update in Sports Medicine: Proc. 2nd Scand. Conf. Sports Med. (1987), p. 16. [20] R.A. Fergurson, P. Krustrup, M. Kajaer, M. Mohr, D. Ball, J. Bangsbo, J. Appl. Physiol. 101 (1) (2006) 47. [21] J. Bangsbo, P. Krustrup, J. Gonzalez-Alonso, B. Saltin, Am. J. Physiol. Endocrinal. Metab. 280 (2001) E956. [22] P. Krustrup, J. Gonzales-Alonso, B. Quistorff, J. Bangsbo, J. Physiol. 536 (2001) 947. [23] P. Krustrup, R.A. Ferguson, M. Kajaer, J. Bangsbo, J. Physiol. 549 (2003) 255. [24] J. Gonzalez-Alonso, B. Quistorff, P. Krustrup, J. Bangsbo, B. Saltin, J. Physiol. 524 (2000) 603. [25] B. Saltin, A.P. Gagge, U. Bergh, J.A. Stolwijk, J. Appl. Physiol. 32 (5) (1972) 635. [26] B. Saltin, in: A. Taylor, P.D. Gollmack (Eds.), Biochemistry of Exercise VII, Human Kinetics, Champaign IL, 1990, p. 387. [27] D. Linnarsson, Acta Physiol. Scand. Suppl. 415 (1974) 1. [28] P.E. Di Prampero, L. Peeters, R. Margaria, J. Appl. Physiol. 34 (1973) 628. [29] J.R. Stirling, M.S. Zakynthinaki, V.L. Billat, Bul. Math. Biol. 70 (5) (2008) 1348. [30] M.S. Zakynthinaki, J.R. Stirling, M. Sillero, J. Sampedro, I. Refoyo, Materials Matemàtics, UAB, 2007. [31] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, second ed., Cambridge University Press, 1993.
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Stochastic optimization for the calculation of the time dependency of the physiological demand during exercise and recovery M.S. Zakynthinaki a,b,∗ , J.R. Stirling b a b
Instituto de Ciencias Matemáticas, CSIC - UAM - UC3M - UCM, c/Serrano 121, 28006 Madrid, Spain Facultad de Ciencias de la Actividad Física y del Deporte, Universidad Politécnica de Madrid, Avd. Martin Fierro s/n, 28040 Madrid, Spain
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 25 January 2008 Received in revised form 3 June 2008 Accepted 31 July 2008 Available online 5 August 2008 PACS: 05.45.Tp 07.05.Tp 47.10.Fg 47.63.Cb 87.19.Hh
The stochastic optimization method ALOPEX IV is successfully applied to the problem of estimating the time dependency of the physiological demand in response to exercise. This is a fundamental and unsolved problem in the area of exercise physiology, where the lack of appropriate tools and techniques forces the assumption and the use of a constant demand during exercise. By the use of an appropriate partition of the physiological time series and by means of stochastic optimization, the time dependency of the physiological demand during heavy intensity exercise and its subsequent recovery is, for the first time, revealed. © 2008 Elsevier B.V. All rights reserved.
Keywords: ALOPEX stochastic optimization Demand Modelling Physiological time series Heart rate Nonlinear dynamical systems
1. Introduction 1.1. A brief introduction to ALOPEX IV The version IV of the stochastic optimization method ALOPEX (ALgorithm Of Pattern EXtraction) has been proven to be the fastest and easiest in its implementation [1–3] when compared to other, previous versions of the same method (versions ALOPEX I, ALOPEX II [3–5] and ALOPEX III [3,6,7]). All ALOPEX methods, however, are shown to be very powerful, especially in optimization problems of many variables in real time. They are fast, effective and very easy to implement. The main advantage of ALOPEX methods is that no knowledge of the dynamics of the system or of the functional dependence of the cost function on the control variables, is required. To give a brief introduction to ALOPEX IV stochastic optimization, let us assume that the aim is to find the values of the N control variables x1 , x2 , . . . , x N that maximize function
f (x1 , x2 , . . . , x N ) = f (x) (the cost function). The method of ALOPEX IV works as follows: if xi (k) is the value of the ith control variable after the kth iteration and f (k) (x1 (k) , x2 (k) , . . . , x N (k) ) is the value of the control function after the kth iteration, then, the value of the ith control variable in the next (k + 1)th iteration is (k+1)
xi
(k)
= xi + c xi (k)
f (k) (k) + gi , | f (k−1) |
i = 1, . . . , N
(1)
where
xi ≡ xi (k) − xi (k−1) , (k)
f (k) ≡ f (k) − f (k−1) , and
f (k−1) ≡ f (k−1) − f (k−2) .
*
Corresponding author at: Instituto de Ciencias Matemáticas, CSIC - UAM - UC3M - UCM, c/Serrano 121, 28006 Madrid, Spain. Tel.: +34 626575554. E-mail addresses: [email protected] (M.S. Zakynthinaki), [email protected] (J.R. Stirling). 0010-4655/$ – see front matter doi:10.1016/j.cpc.2008.07.012
©
2008 Elsevier B.V. All rights reserved.
A proof of convergence of ALOPEX IV, in the absence of noise, can be found in [8]. For more detailed discussion on ALOPEX IV see [3,8].
M.S. Zakynthinaki, J.R. Stirling / Computer Physics Communications 179 (2008) 888–894
1.2. Physiological background Let us denote a physiological variable as s( v , t ) (referring to either the heart rate, or the rate of change of Oxygen uptake), assuming that the current value of s depends on the intensity of the particular exercise (velocity v) and the time passed after the beginning of the exercise, t. The understanding of s( v , t ) in response to exercise, is a fundamental area of exercise physiology and has many applications in sport, see [9–13] as well as medicine and health in general [14]. The values of s( v , t ) always lie within the physiological limits of that particular variable, i.e. smin s( v , t ) smax (for example there is a maximum rate that the heart is able to sustain, as well as a minimum rate, which corresponds to absolute resting values). In the present study we focus on exercises of constant intensity (v constant), as this simplest case is currently the main area of research in exercise physiology. The term ‘on-transient’ kinetics refers to the change in the values of s( v , t ) in response to an increase in the intensity of the exercise, while the term ‘off-transient’ refers to the kinetics of s( v , t ) as the body recovers to a new lower value of s, following a decrease in the exercise intensity. It has been observed that the on transient kinetics of s( v , t ) are slightly different for different constant exercise intensities. For moderate intensities, s( v , t ) follows an exponential-like rise [15,16] until it reaches a steady state. For very heavy intensity exercise the final steady state value is the maximum value smax ; there is, however, a delay (a slowing down of the kinetics) in reaching this value. For severe intensity exercise the values of s( v , t ) rise very steeply and approximately exponentially, until they are limited by the maximum value smax , which then becomes the steady state for the remaining time the exercise can be carried out at. Note that if the exercise is too severe then the subject may have to stop due to fatigue before they reach the smax . Let us assume the velocity and time-dependent function D ( v , t ) 0 that describes the physiological demand for that particular exercise (see also [8,17,18]). It is generally assumed that, for constant intensity exercises where the demand is such that D ( v , t ) smax , the value of the demand is constant and equal to the steady state value that s( v , t ) finally reaches [19] (the asymptote of the physiological time series). This time-independent value D ( v ) depends on the intensity of the particular exercise can easily be obtained from the time series of s( v , t ). For severe or very high intensity exercise there is D ( v , t ) smax . In the literature it is commonly assumed that, even in these cases, the demand is only a function of v and does not depend on time [26]. However this assumption can be shown to be an approximation which is not valid in general [8]. Indeed, for higher values of demand, where D ( v , t ) > smax , the function of demand is probably also a function of time. It is known that the energy cost of accelerating from a particular value of velocity to another value of velocity is greater, due to inertia, than the energy cost of remaining at that velocity. It is also known that the efficiency of muscle decreases with time, during constant intensity exercise [20–24]. For prolonged exercise and at a given constant intensity the energy cost increases due to a loss of efficiency during the exercise. Efficiency has also been shown to be dependent on exercise intensity, with the mechanical efficiency being lower in exhaustive exercise than in sub-maximal exercise [19,25,26]. The on-transient demand should be therefore considered to be a function of time, both in the initial and also the final stages of a bout of heavy intensity exercise. In the case of recovery (off-transient), when the intensity of the exercise is reduced to a new constant level, there is also a time dependency in the demand. It has been observed that in the initial stages of recovery following a very heavy work load, high values of s( v , t ) [27] persist for some seconds before they rapidly drop (see
889
also [28]). The reduction in the heart rate then begins to decrease as we approach the new recovery demand. There is then an ultraslow reduction in the heart rate depending on the severity and duration of the previous exercise; this can last for several hours [9]. It is known that the time to reach resting level of heart rate or oxygen uptake can be very long following very intense exercise, in fact it may not reach resting levels the same day of the exercise [9, 26]. These phenomena suggest that the demand remains elevated for a period of time during recovery. 1.3. The present study A previous study by the authors [8] showed that the ontransient heart rate demand (and also the oxygen demand [29]) can be modelled, using the model of Stirling et al. [17,18], as a constant, only for sufficiently low exercise intensities and for the off-transients following sufficiently low on-transient exercise intensities. The same study showed, however, that these assumptions break down for higher exercise intensities; a model curve assuming constant, time-independent demand could not be optimally fit to the time series of s( v , t ) for heavy intensity exercise [8]. This result justified the expectation for the demand stated in the paragraphs above. The demand, therefore, during exercise and recovery for physiological variables such as the heart rate and oxygen uptake is, in general, a function of time and not a constant, as is usually assumed. The present study applies ALOPEX IV stochastic optimization to the problem of optimizing the fit of the model [17, 18] to physiological time series data (see also [8]) with the aim of revealing, through the process of optimization, the time dependent nature of the physiological demand during exercise and recovery. For the purposes of our study heart rate time series data corresponding to heavy intensity exercise and its subsequent recovery is used. The data consists of two sets of exercise and recovery for two different heavy exercise intensities. This data is the same data that was used in [8] to show that the assumption of a constant demand during heavy intensity exercise and the recovery that follows is incorrect. Our method assumes a partition of the physiological time series in respect to time and calculates a value of the demand for each subset of the data (see Section 3). 2. The optimization problem 2.1. Modelling the kinetics of s( v , t ) The kinetics of a physiological variable s( v , t ) have been recently modelled [8,17,18] using a set of coupled ordinary differential equations of the form
˜ s( v , t ) − smin s˙ ( v , t ) = A v˙ = I (t )
B
C
smax − s( v , t )
E
D ( v , t ) − s( v , t )
(2) (3)
˜ , B , C and E where v is the velocity of the particular exercise, A are parameters that characterize the subject’s current fitness condition (for a discussion on the meaning of these parameters see [8, 17,18]), D ( v , t ) is the demand of the particular exercise (assumed to be dependent on the velocity of the exercise and time) and smin s( v , t ) smax . The function I (t ) describes the rate of change of velocity; in the present study we will assume that I (t ) = 0 and therefore v = constant for a particular exercise. For numerical reasons and without loss of generality, we assume the dimensionless normalized variable s˜ ( v , t ) ≡
s( v , t ) − smin smax − smin
890
M.S. Zakynthinaki, J.R. Stirling / Computer Physics Communications 179 (2008) 888–894
so that 0 s˜ ( v , t ) 1, and a function of a normalized demand
˜ (v , t) ≡ D
that, in accordance with the analysis of the model presented in [18], during the modelling and optimization process of the present work:
D ( v , t ) − smin smax − smin
see also [8,29]. The dynamics of the dimensionless normalized variable s˜ ( v , t ) will then be described by the equation
B
s˙˜ ( v , t ) = A s˜ ( v , t )
C
1 − s˜ ( v , t )
E
˜ ( v , t ) − s˜ ( v , t ) D
(4)
˜ (smax − smin ) B +C + E −1 . where A is the dimensionless parameter A = A Assuming that the parameters A , B , C and E of Eq. (4) are ˜ ( v , t ) that describes the demand, known and so is the function D Eq. (4) can be integrated in respect to time to give the time series {˜sm } N −1 that model the response of the physiological variable k k=0 s˜ (normalized heart rate), see also [8,18,29]. 2.2. Fitting the model to physiological time series data—The optimization parameters
• the value of the parameter E is equal to E = 1 • the value of the parameter C lies within the limits cos(C π ) > 0. Taking the above into account, it is easy to see that the shape of }N −1 that models the data depends on the values of the curve {sm k k=0 the parameters A , B and C of Eq. (4). As mentioned before these values are characteristic of the particular subject, do not change with exercise intensity and give an indication of the subject’s level of fitness, see also [8,17,18,29]. When the values of A , B and C are kept constant for a particular subject, what is left to control the shape of the model curve are the values of D ( v , t ) which, in the present study, are assumed to be time-dependent. For an exercise of a particular constant intensity the control variables of the optimization are the values of the function D (t ).
In the following sections we drop the tildes for simplification. This way the variable s( v , t ) will refer, from now on, to the normalized variable s˜ ( v , t ) and the function D ( v , t ) will refer to the ˜ ( v , t ). normalized demand D
3. The method
Let us assume that the time series describe the basic response pattern of the data (the smooth curve describing the kinetics of the data [30]). For each point k = 0, . . . , N − 1 the sum of the vertical distances (residuals R) between the smooth curve describing the data and the curve {sm } N −1 provided by the model is k k=0
As the aim of the present work is to emphasize the need of time dependent demand, we will use the same data sets used in our previous study [8] that correspond to heart rate data. We note that the experiment was carried out on a tartan track and followed to as good an approximation as possible a square wave protocol consisting of five work periods, of four laps each (for more details on the experimental protocol see [8]). The velocities of the on on sets were, v on 1 = 13.4 Km/h, v 2 = 14.4 Km/h, v 3 = 15.7 Km/h, on v on = 17 . 0 Km / h and v = 17 . 9 Km / h, where the superscript on 4 5 refers to the on-transient exercises. Every exercise was followed by a static recovery of 10 minutes, during which the subject lay hor-
f N −1 {sk }k=0
R=
N −1
f
sk − sm k
2
.
k=0
As described in [8,17], the curve that the model calculates is considered to be an optimal fit to the smoothed data curve {sk }kN=−01 if the value of R is minimum. This is an optimization problem equivalent to the problem [8] of maximizing a cost function f defined as f
f ≡ −R = −
N −1
f
sk − sm k
2
.
(5)
k=0
The aim of the optimally fitting the model curve to the data is achieved via correctly estimating the demand that corresponds to the particular data set and appropriately modifying the shape of the curve provided by the model, i.e. by changing the values of the parameters A , B , C and E of the model [17,18]. We note
3.1. The data
off
izontally and still on the floor (v i = 0, for i = 1, . . . , 5 where the subscript off refers to the off-transient recovery period). In what follows we focus on the study of the heart rate data of exercises 4 and 5, as this data corresponds to heavy intensity exercise. Indicatively, Fig. 1 shows the time series data (raw and normalized) of exercise 4 for the total time of the exercise (338 s on-transient and 600 s off-transient). In this figure, the heart rate data of 20 s before the start of the on-transient part is also shown. Fig. 2 presents, as an indication of the speed time series data, the actual time series of the velocity during the on-transient period of exercise 4. As can be seen, the velocity is as close to
Fig. 1. Raw normalized data, exercise 4. On-transient: 321 s, off-transient: 600 s.
M.S. Zakynthinaki, J.R. Stirling / Computer Physics Communications 179 (2008) 888–894
891
Fig. 2. Exercise 4, on-transient. Velocity time series.
constant as possible in an experimental situation, including measurement error.
the length of the subsets will be too large to consider a constant demand • an M too large (corresponding to a large number of small subsets) will result in miscalculation of D i (tk ): the length of the j subsets will not be enough to correctly calculate a value of D i .
3.2. The cost function The model of Eq. (4) was optimally fit to the data sets described in Section 3.1 above, as described in our previous study [8]. The assumption of a constant, time-independent demand resulted in an optimal fit of the model only for the data sets 1–3, the ones that corresponded to medium exercise intensity. The optimal parameter values, for the particular subject were found to lie within a narrow parameter neighborhood. More specifically, the best parameter combination that optimally fit the model to the data was found to be A = 0.54,
B = 1.63,
C = 1.75.
With these parameter values Eq. (4) becomes
1.63
s˙ ( v , t ) = 0.54 s( v , t )
1.75
1 − s( v , t )
D ( v , t ) − s( v , t ) .
(6)
Let us consider the ith filtered data set [30], where i = 4, 5 as we consider only the data of exercises 4 and 5. This data consists of the N points {si k }kN=−01 and the time takes the discrete values tk , k = 0, . . . , N − 1. The differential, in respect to time, Eq. (6) has the demand D i (tk ) as a variable: f
f 1.63
si˙m = 0.54 si k k
f 1.75
1 − si k
f
D i (tk ) − si k
k = 0, . . . , N − 1, i = 4, 5.
(7)
{s i m } N −1 k k=0
Time integration of (7) yields the time series that provide the model curve. We note here that for the numerical integration, the fifth-order Cash–Karp Runge–Kutta method [31] was used. The time series of the demand D i (tk ) are calculated by the maximization of the cost function as described by Eq. (5). 3.3. The partition of the data sets Let us consider a partition of the data set into M parts. The partition is done in respect to time, such that the time duration (t 1 −t 0 ) . For the jth subset of the data, where of each subset is N −M
By trial and error it was found that the optimal partition was such that the time duration of each data subset was within the interval [18 s, 40 s]. 4. Results 4.1. Estimating the values of the demand from the time series data In the sections that follow we indicatively present the results of a time partition of 20 s. We concentrate on three main areas of the heart rate data: the first part of the on-transient, the very last part of the on-transient and the initial stages of the off-transient and finally, the rest of off-transient time series. As mentioned in Section 3.1, Fig. 1 indicatively shows the time series data of exercise 4, both for the on- and the off-transient. In Fig. 1 we have divided the data into 5 sections: (1) section A includes the first 200 s (approximately) of the ontransient (2) section B includes the second part of the on-transient, from approximately 200 s up to 20 s before the end of the ontransient (3) section C includes the last 20 s of the on-transient (4) section D includes the first 20 s of the off-transient (5) section E includes the remaining of the off-transient data. In Fig. 1 the heart rate data for t ∈ [−20, 0] before the start of the exercise are also shown. It should be noted here that the on-transient part of both exercises 4 and 5 was preceded by a rest period for which it was possible to estimate a constant value of the demand. This value was equal to the asymptotic value of s during rest, which was found to be equal to 0.227 for both exercises 4 and 5.
j
j = 1, . . . , M, we consider a demand of constant value D i that will be calculated via the optimization process described in Section 3.2. The time series D i (tk ) of the demand of the ith data set will then j be composed of all the D i values, j = 1, . . . , M. It should be noted here that the process of finding the demand D i (tk ) is strongly dependent on the choice of M, as
• an M too small (corresponding to a small number of large subsets) will result in loss of detail in the time series of D i (tk ):
4.1.1. Estimating the demand during the first part of the on-transient (area A, Fig. 1) Eq. (8) gives the estimated values of the demand, for the time intervals of the on-transient part of exercise 4,
j on D4
=
⎧ 1.0723, ⎪ ⎪ ⎨ 0.8951, ⎪ 0.9132, ⎪ ⎩ 0.9481,
j=1 j=2 j=3 j4
time partition = 20 s
(8)
892
M.S. Zakynthinaki, J.R. Stirling / Computer Physics Communications 179 (2008) 888–894
and Eq. (9) shows the results for the on-transient part of exercise 5:
j on
D5
=
⎧ 1.1102, ⎪ ⎪ ⎨ 0.9362, ⎪ 0.9435, ⎪ ⎩ 0.9591,
j=1 j=2 j = 3, 4
time partition = 20 s.
(9)
j5
The results were as expected: as can be seen by observation of Eqs. (8) and (9), the demand reaches an highly elevated value straight after the beginning of the exercise. This value drops after about 20 s after the beginning of the exercise and then rises again. 4.1.2. Estimating the demand during recovery (area E, Fig. 1) The estimated values of the demand calculated for the offtransient part of exercise 4 and for a time partition of 20 s are given in Eq. (10) below. It should be noted here that, in Eq. (10), there is j = 1, . . . , 30, since the total time of the off-transient part of each exercise was 600 s. It should also be noted that the values of the demand for the first 20 s of the off-transient are found by a combination of the results of the present section and the results of the following Section 4.1.3.
⎧ 0.2901, j = 1, 2 ⎪ ⎪ ⎪ ⎪ . 2852 , j=3 0 ⎪ ⎪ ⎪ ⎪ ⎪ j = 4, 5, 6 0.2651, j off ⎨ D4 = 0.2453, j = 7, 8, 9, 10 ⎪ ⎪ ⎪ . 2335 , j = 11, 12, 13, 14, 15 0 ⎪ ⎪ ⎪ ⎪ ⎪ 0.2208, j = 16, 17, 18, 19 ⎪ ⎩ j 20 0.2134, time partition = 20 s.
(10)
Eq. (11) gives the estimated values of the demand as calculated for the off-transient part of exercise 5:
j off D5
=
⎧ 0.4294, ⎪ ⎪ ⎪ ⎪ 0.3652, ⎪ ⎪ ⎪ ⎪ 0.3064, ⎪ ⎪ ⎪ ⎨ 0.2827,
j=1 j=2 j=3
j = 4, 5 ⎪ . 2645 , j = 6, 7 0 ⎪ ⎪ ⎪ ⎪ . 2538 , j=8 0 ⎪ ⎪ ⎪ ⎪ ⎪ ≈ 0.2500, 9 j 22 ⎪ ⎩ ≈ 0.2450, j 23
time partition = 20 s.
(11)
As was expected (see Section 1.2) the heart rate demand is found to be elevated at the beginning of the off-transient period and slowly decreasing until it reaches an almost steady state (if the recovery period could last for longer, the values of the demand would asymptotically approach the minimum rest value of D = 0). 4.1.3. Estimating the demand during the last part of the on-transient and the initial stages of the off-transient (areas C and D, Fig. 1) In cases where the exercise intensity is heavy (as is the case for exercises 4 and 5) an estimation of the values of the demand for the last part of the exercise is not possible. It could be said that the last part of the on-transient time series data, where the values of s( v , t ) reach an asymptote (which is near smax , as the intensity is heavy) is a “blind part” of the time series data. In the data of exercise 4 (see Fig. 1) this “blind part” could be said to start after about 200 s from the start of the exercise. Our method estimates the value of the demand at the very last part of the on-transient time series data by making use of the fact that the kinetics of s( v , t ) can be described as a dynamical system, see Eq. (4). When time is allowed to go backwards this system will be therefore governed by the same dynamics. The kinetics backwards in time of the first part of the off-transient time
series will show what demand preceded the off-transient and will this way reveal the heart rate demand at the very last part of the on-transient exercise. We followed the same procedure described in Sections 4.1.1 and 4.1.2; the only difference is that we allowed time to go backwards, i.e. the procedure started from the off-transient and finished at the on-transient. Using the off-transient data of the first 20 s of the off-transient part (area D, Fig. 1) a value for the demand was estimated for the last 20 s of the on-transient time series data (area C , Fig. 1). The estimated (normalized) values of the demand for the last part of the on-transient of exercises 4 and 5 were found to be equal to 1.04 and 1.09 respectively. Note that these values are also the values of the demand for the beginning of the off-transient: during the first 20 s of the off-transient, the demand drops from 1.04 to 0.2901 and from 1.09 to 0.4294 for exercises 4 and 5 respectively (see Eqs. (10) and (11)). 4.2. The time dependency of the demand Figs. 3 and 4 summarize the results presented in the previous Sections 4.1.2, 4.1.1 and 4.1.3. The dashed vertical lines in these figures denote the end of the on-transient and the beginning of the recovery (off-transient) period. A smooth curve has been fitted through the estimated values of the demand, both for the on- and the off-transient part of the exercises. As discussed before, see Section 4.1.3, no information is provided regarding the “blind part” of the time series data (area B, Fig. 1). It is assumed, however, that the values of the demand throughout the exercise follow a smooth increase, as the exercise intensity was constant. In Figs. 3 and 4 we have estimated the time dependency of the demand for this “blind” area by smoothly connecting the estimated values of the demand for the first t 200 s part of the on-transient (area A) with the estimated elevated value of the very last part of the on-transient (area C ), see the dashed curves in Figs. 3 and 4. Indicatively we present in Fig. 5 the (excellent) fit of the model curve on the on-transient time series data of exercise 4. The values of the demand used for the model curve were the ones given in Eq. (8). We also present in Fig. 6 the fit of the model to the off-transient time series of exercise 5 using the values of the demand of Eq. (11). For comparison, Fig. 7 presents the fit of the model when a demand of constant value ( D 5 )off = 0.25 is considered throughout the off-transient period. 5. Conclusions We present here a novel method for estimating the time dependency of the demand during exercise and recovery, which is based on stochastic optimization. In particular, the stochastic optimization method ALOPEX IV was successfully applied to the problem of optimally estimating the values of the demand of a set of heart rate data. The cost function of the optimization was a measure of the optimal fit of a recently developed dynamical systems model [18] on the time series of the data. For this reason, the data was partitioned in respect to time and the value of the demand, for each time interval, was treated as an optimization variable. The same data set was used in a previous study by the authors [8] for the calculation of the optimal values of the model’s parameters. For the study presented in [8] a constant value for the demand was assumed during exercises of constant intensity. The results, however, showed that such an assumption is not valid for exercises of heavy intensity. The study therefore concluded that a time dependent demand should be considered for optimally modelling the kinetics of the physiological time series data (in this case, heart rate data). As can be observed in Figs. 5 and 6 (and
M.S. Zakynthinaki, J.R. Stirling / Computer Physics Communications 179 (2008) 888–894
Fig. 3. Exercise 4. Estimated normalized demand. On-transient: 338 s, off-transient: 600 s.
Fig. 4. Exercise 5. Estimated normalized demand. On-transient: 321 s, off-transient: 600 s.
Fig. 5. Exercise 4. 338 s on-transient heart rate data (normalized) and model curve. Time dependent demand.
893
894
M.S. Zakynthinaki, J.R. Stirling / Computer Physics Communications 179 (2008) 888–894
Fig. 6. Exercise 5. 600 s off-transient heart rate data (normalized) and model curve. Time dependent demand.
Fig. 7. Exercise 5. 600 s off-transient heart rate data (normalized) and model curve. Constant demand.
especially in Fig. 7 where a constant value of the demand is used to model the off-transient kinetics of the exercise), the success of the method we present here is obvious. The novelty of our method is summarized in Figs. 3 and 4 that show the time dependency of the demand for the data of exercises 4 and 5 respectively. The current literature in exercise physiology assumes a constant demand for an exercise of constant intensity, even though there are many errors in this assumption. The present work shows however that, by the use of a simple method such the one we propose here, it is possible to reveal the values of the demand that are hidden in the kinetics of the data. This way, no assumptions on the values of the demand are necessary and therefore new physiological observations or solutions to classical problems in exercise physiology are now possible. Acknowledgements This work was supported by the programs Ramón y Cajal 2004 and I3 2006, Ministerio de Educación y Ciencia, Spain. References [1] M.S. Zakynthinaki, Y.G. Saridakis, Numer. Algorithms 33 (1) (2003) 509. [2] M.S. Zakynthinaki, Y.G. Saridakis, Comp. Phys. Commun. 150 (3) (2003) 274. [3] M.S. Zakynthinaki, Stochastic optimization for adaptive correction of atmospheric distortion in astronomical observation, PhD thesis, Technical University of Crete, 2001. [4] E. Harth, E. Tzanakou, Vision Research 14 (1974) 1475. [5] T. Tzanakou, R. Michalak, E. Harth, Biol. Cybern. 35 (1979) 161. [6] T. Kalogeropoulos, Y.G. Saridakis, M.S. Zakynthinaki, Comp. Phys. Commun. 99 (1997) 255. [7] Y.G. Saridakis, M.S. Zakynthinaki, T. Kalogeropoulos, Internat. J. Appl. Sci. Comp. 5 (3) (1999) 252. [8] M.S. Zakynthinaki, J.R. Stirling, Comp. Phys. Comm. 176 (2) (2007) 98. [9] P.O. Astrand, K. Rodahl, H.A. Dahl, S.B. Stromme, Textbook of Work Physiology: Physiological of Bases of Exercise, Human Kinetics, 2003.
[10] G.A. Brooks, T.D. Fahey, T.P. White, K.M. Baldwin, Exercise Physiology, Human Bioenergetics and its Applications, third ed., McGraw-Hill, 1999. [11] D.E. Martin, P.N. Coe, Better Training for Distance Runners, Human Kinetics, Champaign IL, 1997. [12] T. Noakes, Lore of Running. Discover the Science and Spirit of Running, fourth ed., Human Kinetics, Champaign IL, 2002. [13] D.C. Poole, A. Jones, Oxygen Uptake Kinetics in Sport, Exercise and Medicine: Research and Practical Applications, Routledge, 2004. [14] K. Wasserman, J.E. Hansen, D.Y. Sue, R. Casaburi, B.J. Whipp, Principles of Exercise Testing and Interpretation, Including Pathophysiology and Clinical Applications, third ed., Lippincott Williams and Wilkins, 1999. [15] A.V. Hill, L. Lupton, Q. J. Med. 16 (1923) 135. [16] R. Margaria, F. Mangili, F. Cuttica, P. Cerretelli, Ergonomics 8 (1965) 49. [17] J.R. Stirling, M.S. Zakynthinaki, I. Refoyo, J. Sampedro, J. Nonlin. Math. Phys. (2008), in press. [18] J.R. Stirling, M.S. Zakynthinaki, B. Saltin, Bul. Math. Biol. 67 (5) (2005) 989. [19] B. Saltin, in: S. Maehlum, S. Nilsson, P. Renstrom, (Eds.), An Update in Sports Medicine: Proc. 2nd Scand. Conf. Sports Med. (1987), p. 16. [20] R.A. Fergurson, P. Krustrup, M. Kajaer, M. Mohr, D. Ball, J. Bangsbo, J. Appl. Physiol. 101 (1) (2006) 47. [21] J. Bangsbo, P. Krustrup, J. Gonzalez-Alonso, B. Saltin, Am. J. Physiol. Endocrinal. Metab. 280 (2001) E956. [22] P. Krustrup, J. Gonzales-Alonso, B. Quistorff, J. Bangsbo, J. Physiol. 536 (2001) 947. [23] P. Krustrup, R.A. Ferguson, M. Kajaer, J. Bangsbo, J. Physiol. 549 (2003) 255. [24] J. Gonzalez-Alonso, B. Quistorff, P. Krustrup, J. Bangsbo, B. Saltin, J. Physiol. 524 (2000) 603. [25] B. Saltin, A.P. Gagge, U. Bergh, J.A. Stolwijk, J. Appl. Physiol. 32 (5) (1972) 635. [26] B. Saltin, in: A. Taylor, P.D. Gollmack (Eds.), Biochemistry of Exercise VII, Human Kinetics, Champaign IL, 1990, p. 387. [27] D. Linnarsson, Acta Physiol. Scand. Suppl. 415 (1974) 1. [28] P.E. Di Prampero, L. Peeters, R. Margaria, J. Appl. Physiol. 34 (1973) 628. [29] J.R. Stirling, M.S. Zakynthinaki, V.L. Billat, Bul. Math. Biol. 70 (5) (2008) 1348. [30] M.S. Zakynthinaki, J.R. Stirling, M. Sillero, J. Sampedro, I. Refoyo, Materials Matemàtics, UAB, 2007. [31] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, second ed., Cambridge University Press, 1993.