These authors share last author Brent Baker is the senior scientist from the lab that motivated the problem, and Amy Herring is the chair of the doctoral research committee that guided this research.

The statistics literature on functional data analysis focuses primarily on flexible black-box approaches, which are designed to allow individual curves to have essentially any shape while characterizing variability. Such methods typically cannot incorporate mechanistic information, which is commonly expressed in terms of differential equations. Motivated by studies of muscle activation, we propose a nonparametric Bayesian approach that takes into account mechanistic understanding of muscle physiology. A novel class of hierarchical Gaussian processes is defined that favors curves consistent with differential equations defined on motor, damper, spring systems. A Gibbs sampler is proposed to sample from the posterior distribution and applied to a study of rats exposed to non-injurious muscle activation protocols. Although motivated by muscle force data, a parallel approach can be used to include mechanistic information in broad functional data analysis applications.

^{123}

Studies of physiologic response to muscle stress are important in developing treatment protocols to combat work-, athletic-, and age-related injury. In order to investigate muscle adaptation and maladaptation following repetitive resistance-type exercises, scientists often obtain a series of functional measures (typically at the beginning and end of a multi-session exercise protocol) on the force produced by the muscle as it moves through its range of motion. These force curves can be compared to determine the benefit/harm of an exercise routine to a population of interest.

We investigate one such study conducted in rodents. In this study scientists are interested in differences in physiologic response between young (3 month) and old (30 month) rats exposed to the same resistive exercise protocol. Here rats underwent 13 training sessions, described in

Each observation was recorded as a muscle force tracing, as illustrated in

We have 84 such force tracings, and wish to model the isometric and stretch shortening force generation. The data are defined as follows: for an individual measurement, 565 evenly spaced functional observations were taken. These measurements were taken twice (pre- and post-exercise protocol) resulting in 2 × 565 = 1130 functional measurements per animal. All 42 animals (27 old and 15 young) underwent the same resistive exercise protocol resulting in 47, 460 total measurements.

As the observed function is thought to be the product of two unknown functions (the isometric and stretch shortening components) defined by ordinary differential equations (ODEs), the individual components are not directly identifiable using standard functional approaches that allow each component to be essentially any continuous function. Consequently they must be modeled using simplifying assumptions, but such assumptions may not be interpretable in terms of muscle physiology. Our approach allows these assumptions to based on scientific considerations while retaining the flexibility of the functional approach. Parametric models based upon ODEs do exist but are known to be inadequate for characterizing muscle force tracings (

The force generated by muscle activation, illustrated in

Models for isometric force measurements date back to

In an effort to develop better training/rehabilitation protocols tailored to individual needs, recent studies have investigated how age affects muscle adaptation and maladaptation following specific non-injurious, repetitive, resistance-type loading protocols designed to induce increases in performance and muscle mass. Initial investigations (

Complexities arise when modeling the force tracings of a stretch-shortening contraction. The force output is a product of the isometric force at time

We are interested in modeling

If one assumes muscle force activation is well described by a known ODE there are several recent contributions on Bayesian parameter estimation in this context.

Other parametric work focuses on novel estimation methodologies for non-linear differential equations. For example

From a nonparametric perspective the problem could be approached using a Gaussian process (GP) emulator (

Recent work based upon stochastic differential equations uses differential equations to model data, like the GP emulator approach, but embeds mechanistic information into a stochastic process through stochastic differential equations. In these approaches the forcing function for a differential equation is defined by a stochastic process (the forcing function is the function

Most similar to our approach is the approach of

By examining the reproducing kernel Hilbert space (

Alternatively to finding the analytical solution for the convolution one can use approximation methods to solve the ODE. One approach relies on the MAP-Laplace approximation (

Consider modeling an unknown functional response
_{1}, …, _{n}_{1}, …, _{n}

We incorporate prior information by defining a covariance kernel favoring shapes consistent with mechanistic information specified by differential equations. We assume the information is expressible in the form of a linear ODE
_{0}(_{m−}_{1}(

One can alternatively look at (

In our case, due to the fineness of temporal sampling, the exact solution to the resulting covariance matrix is often extremely ill-conditioned resulting in computational instability. We tried a wide variety of existing methods for addressing ill-conditioning problems in GP regression with no success. The induced covariance of

There is a large literature on approximate solutions to differential equations. Given a set of initial conditions corresponding to _{0}) as well as the first _{0}, Runge-Kutta (RK) methods (see chapter 9 of ^{th}

The Euler-Cauchy approximation recursively defines a solution to _{1}), …, _{n}_{j}_{l}_{−1}). The approximate Euler-Cauchy solution is formed recursively by:
_{l−}_{1}_{l−}_{1}) is a function of the derivative evaluated at _{l−}_{1} and _{l−}_{1} for _{l−}_{1}_{l−}_{1}) = [_{l−}_{1}) + _{0}(_{l−}_{1})_{l−}_{1}] where the function _{0}(·) is a known function described in (_{{}_{l}_{−1,:}} is the previous row. One then defines row

Though we describe the method using the Euler-Cauchy approximation (an ^{2}) approximation method), a similar

One should also investigate the appropriate step size Δ, which also controls the accuracy of the approximation. Too large a step size may result in unacceptable approximation error, and too small step sizes result in unnecessary computational overhead. In our example, with the Euler-Cauchy approximation, our step size was Δ = (1/260). Under these conditions numerical experiments produced results that had a maximum absolute difference on the order of approximately 10^{−3} between the actual and numerical solution across the entire curve.

For the above approximation, sampling from the mechanistic GP proceeds using a series of conditionally conjugate Gibbs steps. The discussion assumes model (^{−1}_{1}(_{1}), …, _{n}_{n}_{1}), …, _{n}_{0}, _{0}). From the above discussion we define _{1}, …, _{z}_{1}(_{m}_{−1}(_{0}, _{0}), which are independent of

Sample r* ~ ^{−1})^{−1} and ^{−1}_{0}, Σ), which is a (_{0} and Σ. Here _{0} is the prior variance of

Sample _{0} + _{0} + (

Marginalizing out r*, i.e., ^{−1}

Gibbs sampling proceeds by iterating over the above steps. To keep the algorithm general, we do not provide an explicit description of sampling of the GP hyperparameters in the above algorithm. However, we have observed the best mixing when these hyperparameters are block sampled with

The mechanistic GP is not directly applicable to the muscle force application, which has the additional complication of decomposing

We define an ODE for

Placing a linearity assumption on _{0}(_{s}_{s}_{s}_{′} be a linear spline on the interval, 0 prior to, and 1 after the SS contraction, which defines a linear interpolation of the function on this range. In order to model a flexible curve we use a large number of evenly spaced splines to estimate

When the joint is moved through its range of motion, the force on the joint is related to the angle of the joint and other factors. Angular motion is often described using a second order differential equation, and we follow this approach. As the exact form of the differential equation is unknown (damping constant) we model this function through the fully specified second order differential equation:

It is further known that the multiplicative effect of _{a}_{b}_{a}_{b}

The RK approximation is used to sample both _{0})′ _{0}, _{0}) and _{0}_{1})′ ~ _{1}, _{1}), initial value vectors for _{1}), …, _{n}_{1}), …, _{n}

In sampling _{1}, …, _{n}_{0}_{1}, …, _{S}_{β}

Putting the prior _{0}, _{0}) over the initial conditions, define _{0}, Σ_{β}^{−1})^{−1} and ^{−1}

We also modify algorithm 1 to sample g* given _{a}_{b}_{a}_{0} and _{b}_{{n,:}}g*, where _{{n,:}} is the last row of

Calculating _{a}_{b}

On the interval [_{a}, t_{b}_{l}

For each _{ijk}

For each

For

For each _{j}

Sample

Inference on the group average curves

We extend the hierarchical mechanistic process to our application. Here
_{ijk}_{ijk}_{ijk}_{ijk}_{ijk}_{ijk}

For _{ijk}_{ijk}_{a}_{ijk}_{b}

Hierarchical extensions in modeling _{ijk}

We conduct a simulation experiment based upon the model developed in (_{ijk}_{ijk}_{ijk}_{ijk}_{j}

We chose weakly informative priors for all hyper parameters. We place a GP prior over _{ijk}^{2}^{2}). We set ^{−2} ~ _{ijk}_{ijk}_{j}_{j}_{ijk}

We collected 50, 000 MCMC samples disregarding the first 5, 000 as a burn-in. Examinations of trace plots for the quantities of interest, the individual curves and curves in the hierarchy, showed adequate mixing. Geweke’s diagnostic tests on the chains (

For the quantities of interest (i.e., the individual level functions and the first level of the hierarchy _{ijk}_{ijk}

With the goal of investigating the effect of non-injurious, repetitive muscle contractions on muscle force generation, we apply our approach to data compiled from

_{ijk}

The model also allows one to look at individual estimates between curves. Though the animals showed no significant differences at the group level for isometric force generation, individual differences were seen.

Note that there is an additional advantage of modeling the latent forcing function as it may be used to generate or support hypotheses. For example, for _{ijk}_{ijk}

This article proposes a flexible nonparametric Bayesian method that takes into account prior scientific information based upon an ODE. This method relies on a sampling algorithm using a Runge-Kutta approximation to the ODE. The nature of the approximation allows for a hierarchical specification at the population level estimates by modeling the latent forcing function directly. The goal was to develop an inferential framework for muscle force tracings, and investigate the effect of non-injurious resistive exercise protocols on muscles of different ages (i.e., young and old muscles). It was important to accurately model both the overall functional response and the two constituent functions, which themselves have scientific interest.

For our specific problem the model is able to show that the dynamic force generated through the stretch shortening contraction may be more informative and specific in showing adaptation and maladaptation following non-injurious mechanical loading. Both groups of rats have an adaptive response in dynamic muscle force produced. Force generated is, at the population level, greater after the exercise protocol, but both groups have statistically non-significant responses in isometric force generation. As maximal isometric force generation, which is represented by a single value, has been seen as the gold standard response of interest in the muscle force literature in terms of measuring adaptation/maladaptation. The models presented here may allow researchers to begin to model more sensitive endpoints and look at their data from a functional data perspective. Further, such analyses may lead to new insights. For example, the maximal force for older rats appears to occur at a different joint angle, suggesting impacts affecting the length=tension curve dynamics as the animals age.

The proposed approach can be extended to human muscle force tracings, and may allow for in-depth study of human physiologic responses to exercise routine post training. Such an analysis, on the entire force output, has previously not been attempted. For such a study, age, as well as other variables, can be included in the hierarchical framework. In this manner the efficacy of the current ‘one size fits all’ approach across the spectrum of prevention/intervention including occupational, athletic, physical therapy, strength and conditioning, and wellness programs can be studied. If similar results are seen in human populations, this might suggest that different treatment protocols are appropriate depending on age or other variables of interest.

In domains outside of the muscle force literature, the approximation method allows researchers to easily incorporate mechanistic information into analyses, while bypassing the cumbersome integrals that are required for direct calculation of the covariance kernel. In situations where there are fewer observations that are not evenly spaced, one can still use this approach. This can be done either by using our regression based approach directly, or, when inference is desired on less closely spaced observations, one can approximate the covariance matrix by computing

The authors would like to thank John Bailer and Rich Lovering as well as three anonymous referees for their comments and suggestions on an earlier version of this manuscript.

This research was supported in part by NIH grants ES017436 and 5R01ES020619.

The findings and conclusions in this report are those of the authors and do not necessarily represent the views of the National Institute for Occupational Safety and Health.

These authors share last author Brent Baker is the senior scientist from the lab that motivated the problem, and Amy Herring is the chair of the doctoral research committee that guided this research.

Fifteen muscle force measurements (gray lines) made on young rats before training, with the mean of these measurements shown in black. The five sections separated by a vertical line describe different portions of the muscle activation. The first and last sections represent force generation when the muscles are not contracting. The second and fourth sections represent the force generated during an isometric contraction; with the third section denoting the stretch shortening contraction, which takes place during an isometric contraction.

Estimated group level curves for the dynamic force in a stretch shortening contraction for a simulated data set. The black line and corresponding 95% credible region (dotted lines) represent the mean estimated curve and credible region. The true curve is given by the gray line and can be seen to be well represented by the model.

The estimated mean isometric force generated for a single animal pre- and post-treatment. The dark black line represents central estimates of the function

The estimated group level dynamic force multiplier generated by old (left column) and young (right column) animals. Here the top row represents the average force increase during a stretch shortening contraction pre- (black line) and post- (dotted line) exercise protocol. The bottom row represents the observed increase (solid black lines) thought to be due to exercise along with the 95% pointwise credible interval (dotted line) for this increase due to exercise.

Estimated mean isometric muscle force generated for the young animals pre (dashed line) and post (solid line). The bottom row gives the estimates, and 95% pointwise credible intervals of the difference between the two.

Estimated dynamic muscle force for an old animal. This figure shows the central estimate for the pre- (dashed line) and post- (solid line) for the isometric contraction of an individual old animal. Credible intervals are not shown as they are indistinguishable from the central estimate.

This table presents the steps of a resistive-type exercise protocol administered in a rat model of muscle function. Degrees indicate the angle of the joint during application of force. We investigate differences in observation 5 from the first and last experiment. Adapted from

Observation | Experimental Test Protocol (3 Times/wk for 4.5 weeks) |
---|---|

Single Isometric force test conducted at 90° | |

1 Stretch-Shortening Contraction at 70° – 140° – 70° | |

8 sets of 10 Stretch-Shortening Contractions 90° – 140° – 90° | |

Single Isometric test conducted at 90° | |

1 Stretch-Shortening Contraction at 70° |