Free
Original Contribution  |   December 2008
Mathematical Analysis of Applied Loads on Skeletal Muscles During Manual Therapy
Author Notes
  • From the Departments of Biomedical Engineering (Drs Chaudhry and Findley) and Mathematical Sciences (Dr Bukiet) at the New Jersey Institute of Technology in Newark, and the War-Related Illness and Injury Study Center at the Veterans Affairs Medical Center in East Orange, NJ (Drs Chaudhry and Findley). 
  • Address correspondence to Bruce Bukiet, PhD, Department of Mathematical Sciences, Center for Applied Mathematics and Statistics, New Jersey Institute of Technology, Newark, NJ 07102-1982. E-mail: bukiet@m.njit.edu 
Article Information
Neuromusculoskeletal Disorders
Original Contribution   |   December 2008
Mathematical Analysis of Applied Loads on Skeletal Muscles During Manual Therapy
The Journal of the American Osteopathic Association, December 2008, Vol. 108, 680-688. doi:10.7556/jaoa.2008.108.12.680
The Journal of the American Osteopathic Association, December 2008, Vol. 108, 680-688. doi:10.7556/jaoa.2008.108.12.680
Abstract

Context: Because of the lack of accurate values for mechanically applied forces in osteopathic manipulative treatment, osteopathic physicians must apply mechanical forces intuitively. However, excessive loading and high velocity maneuvers carry risks for patients.

Objective: To determine the loads required to produce compression, shear, extension, and twist on biceps muscle during manual therapy.

Methods: A mathematical analysis valid for the in vivo state of biceps muscle was performed to determine the loads produced in a simple elastic biceps muscle model and a more realistic viscoelastic biceps muscle model.

Results: Loads of 7% lesser pressure were needed to produce 10% deformation of biceps muscle using the viscoelastic model, compared with the elastic model. In the viscoelastic model, there was stress relaxation of 18% of maximum pressure when muscle was deformed by 10% over 60 seconds and maintained in that state for 200 seconds. With quick maneuvers, the viscoelasticity effect was decreased.

Conclusion: The biceps muscle is 15 times stiffer in the direction parallel to the muscle fibers than in the direction perpendicular to the fibers. The results of the present study may be used by osteopathic physicians to adjust their manual techniques to match viscoelastic properties of specific tissues. Because biceps muscle is viscoelastic, the results obtained with the viscoelastic model would be more useful than the results obtained with the elastic model for determining viscoelastic loads on this muscle.

Several forms of manual therapies, including various osteopathic manipulative treatment (OMT) techniques, have been developed to improve postural alignment and other expressions of musculoskeletal dynamics in patients.1-3 These techniques involve shear, longitudinal stretch extension, deep pressure compression, and unwinding twisting movements.3 
However, excessive mechanical loading and high velocity maneuvers carry some risks for patients and should be approached with caution by manual practitioners.3 At present—because of a lack of known, accurate values for mechanically applied forces—mechanical forces must be applied intuitively, which may be helpful in some cases but harmful in others. Therefore, there exists a need to determine scientifically accurate values of mechanically applied forces used in manual therapy. These values should be based on scientifically established in vivo mechanical properties of soft tissues, such as skeletal muscles and fasciae, so that any manual therapy technique is beneficial to the patient rather than counterproductive. 
To most effectively apply treatment, it is important that in vivo mechanical properties of human skeletal muscles be known a priori. Recently, magnetic resonance elastography (MRE) has been used to quantify the mechanical properties of soft biologic tissues.4,5 The shear wave patterns observed with this technique can be used noninvasively to evaluate the anisotropic elastic coefficients of skeletal muscles.6-8 Anisotropic materials have different physical properties in different directions. Isotropic materials have the same physical properties in all directions. Transversely isotropic materials have a plane (ie, the plane of isotropy) in which the material properties are the same in different directions, but in the direction perpendicular to this plane (ie, the axis of transverse isotropy), the material properties are different than in the plane (Figure 1). 
Skeletal muscle is considered to be transversely isotropic for the analysis of experimental shear wave speeds.7 Therefore, assuming biceps muscle to be a transversely isotropic material, Papazoglou et al9 quantitatively evaluated the mechanical properties of biceps muscle. In their experiment, the actuator was applied at the distal tendon of the biceps muscle to induce shear vibration along the principal axis of symmetry (ie, parallel to the direction of the muscle fibers). The resulting two-dimensional waveform profiles of the wave images were used to directly determine the shear modulus of the biceps muscle as 54±3 kPa.9 The shear modulus is the ratio of shear stress (ie, the tangential stress on a plane) to shear strain (ie, the angle of distortion caused by the shear stress).9 
Papazoglou et al9 used the linear theory of transverse isotropy with small-deformation strain tensor developed by Lai et al10 to determine the Young's moduli along the axis of symmetry and along axes in the plane perpendicular to the axis of symmetry. The linear theory was also used to determine the shear moduli along these axes. According to the linear theory, stress is linearly proportional to strain, and strain is the change in length divided by the original length.10 Young's modulus along a particular direction is defined as the ratio of the stress (ie, the force per unit area) along that direction to the strain along that same direction.10 
Figure 1.
Three-dimensional model of transversely isotropic biceps muscle. Transversely isotropic materials have a plane (ie, the plane of isotropy) in which the material properties are the same in different directions, but in the direction perpendicular to this plane (ie, the axis of transverse isotropy), the material properties are different than in the plane. In this model of biceps muscle, the X1, X2 plane is the plane of isotropy. The elastic properties along the X3 axis, which is the axis of transverse isotropy, are different from the elastic properties in the X1, X2 plane. (Adapted from figure provided by Ingolf Sack, PhD, of the Institute of Radiology at Universitätsmedizin Berlin, Germany, and previously published online at: http://www.elastography.de/.)
Figure 1.
Three-dimensional model of transversely isotropic biceps muscle. Transversely isotropic materials have a plane (ie, the plane of isotropy) in which the material properties are the same in different directions, but in the direction perpendicular to this plane (ie, the axis of transverse isotropy), the material properties are different than in the plane. In this model of biceps muscle, the X1, X2 plane is the plane of isotropy. The elastic properties along the X3 axis, which is the axis of transverse isotropy, are different from the elastic properties in the X1, X2 plane. (Adapted from figure provided by Ingolf Sack, PhD, of the Institute of Radiology at Universitätsmedizin Berlin, Germany, and previously published online at: http://www.elastography.de/.)
In the present study, the quantitative information about the anisotropic elastic properties of biceps muscle obtained by Papazoglou et al9 was used to determine the magnitude of both the elastic and viscoelastic stresses produced in transversely isotropic biceps muscle when subjected to compression, shear, longitudinal extension, and twist in manual therapy. It was hoped that the present analysis would provide useful information to biophysicists for a greater understanding of mechanical properties of muscle, and to osteopathic physicians for beneficial results in patient treatment. 
Methods
To determine the elastic and viscoelastic stresses produced in biceps muscle subjected to compression, shear, extension, and twist, the linear stress-strain relation for transversely isotropic material developed by Lai et al10 as applied to transversely isotropic biceps muscle. This protocol was deemed appropriate because the present study considers small deformations in the muscle. Then, the appropriate mechanical forces needed for producing a specific nondamaging deformation in the biceps muscle were mathematically determined, considering the biceps muscle as both elastic and viscoelastic. 
Deformation
In this section of the study, methods used in evaluating the force (ie, applied load) required to obtain a given pattern of muscle displacement are presented. If an evaluated applied load lies outside the normal physiologic range, then the pattern of displacement corresponding to that load is unlikely to occur during manual therapy. We assumed that the displacements u1, u2, and u3 along the X1, X2, X3 axes were produced on the surface muscle element, covering its entire thickness, by manual therapy techniques involving shear and extension along the X3 axis; compression along the negative X1 axis; and twist about X3 axis (Figure 2). These displacements are given by the following:  
\[\ \begin{array}{rl}u_{1}&=(k_{3}-1)X_{1}+X_{1}(Cos{\psi}-1)-X_{2}Sin{\psi}\\u_{2}&=(k_{2}-1)X_{2}+X_{2}(Cos{\psi}-1)+X_{1}Sin{\psi}\\u_{3}&=k_{1}X_{1}+k_{4}X_{3}\end{array}\]
(1)
 
In this equation, the Xi axes are in the Cartesian coordinate system. In this case, X3 is parallel to the direction of the muscle fibers. The symbol k1 (k1>0) denotes the shear ratio resulting from the application of the tangential force. The maximum shear occurs at the surface of the muscle, where muscle thickness is maximum, and the shear equals 0 at the bottom of the muscle, where there is no muscle thickness. The symbol k2 (k2>0) is the extension ratio resulting from the compression and shear on the surface of the muscle. The symbol k3 (k3<1) denotes the compression ratio resulting from the applied normal pressure. The symbol k4 (k4>0) is the extension ratio resulting from the applied longitudinal force. The symbol ψ is the angle (in radians) caused by the twist applied to the muscle. The derivation of the second terms in the u1 and u2 portions of the equation, corresponding to twist, is provided in the Appendix. 
Figure 2.
Three-dimensional model of skeletal muscle element used in analysis of displacement caused by applied loads. OX1, OX2, and OX3 are the axes in the Cartesian coordinates system; F denotes the longitudinal stress; M denotes the moment; N and T denote the normal pressure and tangential force, respectively. The six faces of this model are designated by combinations of the capital letters A, B, C, D, E, G, H, and I. The symbols a, b, -b, and c represent the coordinate values of the faces along the three axes.
Figure 2.
Three-dimensional model of skeletal muscle element used in analysis of displacement caused by applied loads. OX1, OX2, and OX3 are the axes in the Cartesian coordinates system; F denotes the longitudinal stress; M denotes the moment; N and T denote the normal pressure and tangential force, respectively. The six faces of this model are designated by combinations of the capital letters A, B, C, D, E, G, H, and I. The symbols a, b, -b, and c represent the coordinate values of the faces along the three axes.
The strain, Ei,j, in the linear theory of elasticity is given by the following:  
\[\ E_{\mathrm{i},\mathrm{j}}=\frac{1}{2}(u_{\mathrm{i},\mathrm{j}}+u_{\mathrm{j},\mathrm{i}})\]
(2)
In equation (2), the comma denotes the partial derivative, and i,j=1, 2, 3. 
Using equations (1) and (2), we obtain the following:  
\[\ \begin{array}{rl}E_{11}&=k_{3}-1+(Cos{\psi}-1)\\E_{22}&=k_{2}-1+(Cos{\psi}-1)\\E_{33}&=k_{4}\\E_{12}&=E_{21}=E_{23}=E_{32}=0\\E_{13}&=E_{31}=\frac{k_{1}}{2}\end{array}\]
(3)
Henceforth, we shall designate the term Cosψ - 1 (corresponding to twist) as k5. 
Elastic Skeletal Muscle (Elastic Model)
In this section, the elastic model is used to evaluate elastic stresses that will be used in the equations for the viscoelastic model. The loads required to produce the type of deformation represented by equations (1) through (3) can be determined by evaluating the stresses in the biceps muscle. Stress is defined as the load divided by the area of the cross-section on which the load is applied. 
As previously stated, the biceps muscle is assumed to be a transversely isotropic material (Figure 1).9 The stress-strain relation for transversely isotropic material that is valid for biceps muscle tissue10 is given by the following:  
\[\ \left[\begin{array}{c}T_{11}\\T_{22}\\T_{33}\\T_{23}\\T_{31}\\T_{12}\end{array}\right]=\left[\begin{array}{cccccc}C_{11}&C_{12}&C_{13}&0&0&0\\C_{12}&C_{11}&C_{13}&0&0&0\\C_{13}&C_{13}&C_{33}&0&0&0\\0&0&0&2Y_{13}&0&0\\0&0&0&0&2Y_{13}&0\\0&0&0&0&0&2Y_{12}\end{array}\right]\left[\begin{array}{c}E_{11}\\E_{22}\\E_{33}\\E_{23}\\E_{31}\\E_{12}\end{array}\right]\]
(4)
The following values can be given for the variables in equation (4):  
\[\begin{array}{rl}&C_{11}=\left(\frac{Y_{1}}{1+v_{21}}\right)\frac{(1-v_{31}^{2}(Y_{1}{/}Y_{3}))}{D}\\&C_{12}=\left(\frac{Y_{1}}{1+v_{21}}\right)\frac{(v_{21}+v_{31}^{2}(Y_{1}{/}Y_{3}))}{D}\\&C_{13}=\frac{Y_{1}v_{31}}{D},C_{33}=Y_{3}+2v_{31}C_{13}\\&D=1-v_{21}-2v_{31}^{2}\left(\frac{Y_{1}}{Y_{3}}\right){>}0\end{array}\]
For biologic tissues, we assume the following values9,11:  
\[\begin{array}{rl}&v_{31}=v_{32}=0.45,{\ }\mathrm{for}{\ }D{>}0\\&v_{12}=\left(1-\frac{v_{31}Y_{1}}{Y_{3}}\right)\\&C_{33}=C_{13}\end{array}\]
 
In the above equations, Cij are the elastic coefficients and vij are Poisson's ratios. When a given load is applied along the longitudinal direction, the Poisson's ratio is defined as the ratio of the lateral strain (ie, the strain normal to the longitudinal strain) to the longitudinal strain.12 The symbols Tij and Eij are stress and strain components. The symbols Y1 and Y3 are Young's moduli in the X1 and X3 directions, respectively, and Y12 and Y13 are the shear moduli. The symbol v31 is the Poisson's ratio for the transverse strain in the X1 or X2 direction when stressed in the X3 direction. Similarly, v12=v21 is the Poisson's ratio in the plane of isotropy. 
Thus, from equations (3) and (4), the stresses Tij are given by the following:  
\[\ \begin{array}{rl}T_{11}&=C_{11}E_{11}+C_{12}E_{22}+C_{13}E_{33}\\T_{22}&=C_{12}E_{11}+C_{11}E_{22}+C_{13}E_{33}\\T_{33}&=C_{13}E_{11}+C_{13}E_{22}+C_{13}E_{33}\\T_{23}&=T_{32}=2Y_{13}E_{23}=0\\T_{31}&=T_{13}=2Y_{13}E_{13}\\T_{12}&=T_{21}=2Y_{12}E_{12}=0\end{array}\]
(5)
 
Because all these stresses are constants, equations of equilibrium are automatically satisfied. Next, we consider the following cases: 
  • compression and shear only
  • longitudinal extension only
  • twist and compression only
For the case of compression and shear along the negative X1 axis, we must have T22=T33=0. Then, from equation (5), we get the following:  
\[\ \begin{array}{rl}E_{22}&=\frac{C_{13}-C_{12}}{C_{11}-C_{13}}E_{11}\\E_{33}&=\frac{C_{12}-C_{11}}{C_{11}-C_{13}}\end{array}\]
(6)
The following two equations can then be derived:  
\[\ T_{11}=\left[C_{11}+C_{12}\left(\frac{C_{13}-C_{12}}{C_{11}-C_{13}}\right)+C_{13}\left(\frac{C_{12}-C_{11}}{C_{11}-C_{13}}\right)\right]E_{11}\]
(7)
 
\[\ T_{13}=2Y_{13}E_{13}\]
(8)
 
For the case of longitudinal extension along the X3 axis, we must have T22=T11=0. Then, from equation (5), we derive the following:  
\[\ E_{22}=E_{11}=\frac{-C_{13}E_{33}}{C_{11}+C_{12}}\]
(9)
Equation (9), in turn, leads to the following:  
\[\ T_{33}=C_{13}\left[1-\frac{2C_{13}}{C_{11}+C_{12}}\right]E_{33}\]
(10)
 
For the case of twist about the X3 axis and compression on the faces X1=a, X2=b, we must have T33=0. Then, from equations (3) and (5), we derive the following:  
\[\ \begin{array}{rl}&E_{33}=-(E_{11}+E_{22})\\&E_{11}=E_{22}\\&k_{1}=0\end{array}\]
(11)
 
\[\ \begin{array}{rl}T_{11}&=(C_{11}-C_{13})E_{11}+(C_{12}-C_{13})E_{22}\\T_{22}&=(C_{12}-C_{13})E_{11}+(C_{11}-C_{13})E_{22}\\T_{33}&=T_{13}=T_{23}=0\end{array}\]
(12)
 
The twisting moment, M, applied to the muscle is given by the following equation:  
\[M={{\int}}{{\int}}T_{11}{\cdot}X_{2}dX_{2}dX_{3}+{{\int}}{{\int}}T_{22}dX_{1}dX_{3}\]
The limits of X1 are from 0 to a; of X2, from -b to +b; and of X3, from 0 to c. Thus, M per unit area can be expressed by the following, in which a is the thickness of the muscle:  
\[\ M=\frac{T_{22}a}{2}\]
(13)
 
Viscoelastic Skeletal Muscle (Viscoelastic Model)
The stresses evaluated in the previous section for elastic biceps muscle will now be used to evaluate dynamic stresses for viscoelastic biceps muscle. The elastic model is a simplified model. However, the results derived from the elastic model allow us to evaluate the viscoelastic model, which more closely approximates human tissue. 
To determine the dynamic stress responding to the applied dynamic extension for viscoelastic fasciae, a constitutive equation13 was used to address the viscoelastic behavior in biologic soft tissues. This equation can be written as the following:  
\[\ T_{\mathrm{ij}}(x,t)={{\int}_{-{\infty}}^{\mathrm{t}}}G_{\mathrm{ijkl}}(t-{\tau})\frac{{\partial}T_{\mathrm{kl}}^{(\mathrm{e}){\ }(x,{\tau})}}{{\partial}{\tau}}d{\tau}\]
(14)
In equation (14), Te is the pseudo-elastic stress and a function of strain, which, in turn, is a function of position, x, and time, t; Gijkl (i,j,k,l= 1, 2, 3, 4, 5, 6) is the tensor of relaxation function; τ is a dummy variable for the integration from the beginning (t=0) of the deformation to the present time, t. 
At this stage, we can designate T1, T2... T6 for T11, T22, T33, T23, T31, T12, and E1, E2... E6 for E11, E22, E33, E23, E31, E12. 
Then, equation (14) can be rewritten as the following:  
\[\ T_{\mathrm{n}}={{\int}_{0}^{\mathrm{t}}}G_{\mathrm{nq}}(t-{\tau})\frac{{\partial}T^{(\mathrm{e})}}{{\partial}{\tau}}d{\tau}{\ }(n,q=1,2,3,4,5,6)\]
(15)
In equation (15), the uninvolved diagonal relaxation functions are assumed to be 0, and we retain G11, G22, G33, G13. Then, we get the following:  
\[\ \begin{array}{rl}T_{1}&={{\int}_{0}^{\mathrm{t}}}\left(G_{11}(t-{\tau})\frac{{\partial}T_{1}}{{\partial}{\tau}}+G_{13}(t-{\tau})\frac{{\partial}T_{3}}{{\partial}{\tau}}\right)d{\tau}\\T_{2}&={{\int}_{0}^{\mathrm{t}}}\left(G_{22}(t-{\tau})\frac{{\partial}T_{2}}{{\partial}{\tau}}\right)d{\tau}\\T_{3}&={{\int}_{0}^{\mathrm{t}}}\left(G_{31}(t-{\tau})\frac{{\partial}T_{1}}{{\partial}{\tau}}+G_{33}(t-{\tau})\frac{{\partial}T_{3}}{{\partial}{\tau}}\right)d{\tau}\end{array}\]
(16)
 
The final step in developing field equations for viscoelastic skeletal muscle is to determine the relaxation function. Fung14 introduced a one-dimensional standard model to derive a linear reduced relaxation function with a continuous spectrum. This widely accepted reduced relaxation function14 is presented as the following:  
\[\ G(t)=\frac{1+{{\int}_{0}^{{\infty}}}S(q)e^{-\mathrm{t}{/}\mathrm{q}}dq}{1+{{\int}_{0}^{{\infty}}}S(q)dq}\]
(17)
In equation (17), the relaxation spectrum S(q) denotes the amplitude of the viscous dissipation subject to the frequency 1/q. For many biologic materials that exhibit viscoelastic behavior and low sensitivity to strain rates, the designation of the spectrum with constant amplitude over a range of frequency has been found to fit the experimental results fairly well.14 This spectrum assumes the following form:  
\[\ \begin{array}{rl}S(q)&=\frac{C0}{q},\mathrm{for}{\ }q_{1C}{\leqslant}q{\leqslant}q_{2C}\\S(q)&=0,\mathrm{otherwise}\end{array}\]
(18)
In equation (18), the parameters C0, q1C, and q2C can be determined by experiments examining the stress relaxation of the tissues under constant strain and by experiments of dynamic loading. In the present study, we adopted the parameters obtained from the study of subcutaneous tissues of rats by Iatridis et al.15 Those researchers15 reported that C0=0.25, q1C=1.86 seconds, and q2C=110.4 seconds. Combining equations (15) through (18), we can numerically resolve the viscoelastic stress-time and stress-strain relations for the biceps muscle. 
Because the viscosity is uniform in the muscle tissues, we assume, in equation (18), that G11(t - τ) = G22(t - τ) = G33(t - τ) = G13(t - τ)= G(t - τ), and that the values of the parameters C0=0.25, q1C=1.86 seconds, and q2C=110.4 seconds in G(t - τ) are as in equation (14). 
Results
The results obtained from the above equations will now be used to determine the forces that must be applied, over a set duration, to produce a given amount of stretch in the biceps muscle, assuming that the muscle is either elastic or viscoelastic. 
Because the elastic parameters, such as the Young's moduli and shear modulus for transversely isotropic biceps muscle, have previously been determined by Papazoglou et al,9 we used those parameters to evaluate the mechanical stresses generated by manual therapy on biceps muscle. The values of these parameters are given below:  
\[\ \begin{array}{rl}&Y_{1}=13{\pm}4.5{\ }\mathrm{kPa}\\&Y_{3}=185{\pm}60{\ }\mathrm{kPa}\\&Y_{13}=54{\pm}3{\ }\mathrm{kPa}\end{array}\]
(19)
In equation (19), Y1 and Y3 are Young's moduli along the X1 and X3 directions, respectively (ie, perpendicular and parallel to muscle fibers, respectively). The symbol Y13 is the shear modulus parallel to muscle fibers. 
The following linear and parabolic regime simulations for strain, E, are used to evaluate the forces for the elastic and viscoelastic models: Linear regime:  
\[\ E=\frac{At}{T_{0}}(\mathrm{ie,\ strain\ increases\ linearly\ wiht\ time},t)\]
(20)
Parabolic regime:  
\[\ E=A\left[1-\left(\frac{t-T_{0}}{T_{0}}\right)^{2}\right](\mathrm{ie,\ strain\ varies\ parabolically\ with\ time},t)\]
(21)
In equations (20) and (21), T0 is the time during which treatment is performed in one stroke, and A is the maximum strain, which is attained at the end of that time. 
To evaluate the stresses in equation (16) for the viscoelastic model, we used the stress values as indicated in the parentheses of the integrals in equation (16), which were obtained from the elastic model. Table 1 presents the results of applied forces (ie, combined compression and shear, longitudinal extension alone, combined twist and compression) in terms of stresses for both models, because, in a clinical situation, the area of cross-section and the thickness of the muscle would vary with each patient. 
Table 1
Stress Values for Applied Forces in Elastic and Viscoelastic Models


Stress Value
Applied Force*
Elastic Model
Viscoelastic Model
Compression and Shear, kPa (lb/in2) (k1=0.1; k3=0.9; k5=0)
□ Compressive stress -1.30 (-0.19) -0.80 (-0.12)
□ Shear stress10.80 (1.57)6.64 (0.96)
Longitudinal Extension, kPa (lb/in2) (k1=0; k3=1.0; k4=0.1; k5=0)
□ Longitudinal stress18.50 (2.68)11.40 (1.65)
Twist and Compression, kPa-cm§ (k1=0; k2=0.9; k3=0.9; ψ=10 degrees)
□ Twisting moment23.7014.60
 Abbreviations: k1 denotes shear ratio resulting from application of tangential force; k2, extension ratio resulting from compression and shear on surface of muscle; k3, compression ratio resulting from applied normal pressure; k4, extension ratio resulting from applied longitudinal force; k5, Cosψ - 1; ψ, angle caused by twist applied to muscle.
 *Applied 10% compression and applied 10% shear are on face X1=a of skeletal muscle element. Applied extension is parallel to fiber direction on face X3=c of skeletal muscle element. See Figure 2.
 k2 and k4 can be determined from equation (6) in text.
 k2 can be determined from equation (9) in text.
 §k4 can be determined from equation (11) in text.
 Twisting moment per unit area about the X3 axis.
Table 1
Stress Values for Applied Forces in Elastic and Viscoelastic Models


Stress Value
Applied Force*
Elastic Model
Viscoelastic Model
Compression and Shear, kPa (lb/in2) (k1=0.1; k3=0.9; k5=0)
□ Compressive stress -1.30 (-0.19) -0.80 (-0.12)
□ Shear stress10.80 (1.57)6.64 (0.96)
Longitudinal Extension, kPa (lb/in2) (k1=0; k3=1.0; k4=0.1; k5=0)
□ Longitudinal stress18.50 (2.68)11.40 (1.65)
Twist and Compression, kPa-cm§ (k1=0; k2=0.9; k3=0.9; ψ=10 degrees)
□ Twisting moment23.7014.60
 Abbreviations: k1 denotes shear ratio resulting from application of tangential force; k2, extension ratio resulting from compression and shear on surface of muscle; k3, compression ratio resulting from applied normal pressure; k4, extension ratio resulting from applied longitudinal force; k5, Cosψ - 1; ψ, angle caused by twist applied to muscle.
 *Applied 10% compression and applied 10% shear are on face X1=a of skeletal muscle element. Applied extension is parallel to fiber direction on face X3=c of skeletal muscle element. See Figure 2.
 k2 and k4 can be determined from equation (6) in text.
 k2 can be determined from equation (9) in text.
 §k4 can be determined from equation (11) in text.
 Twisting moment per unit area about the X3 axis.
×
The stress relaxation plots for combined compression and shear of 10% are presented as linear and parabolic strain trajectories in Figure 3. The stress relaxation plots for 10% extension only are given in Figure 4. The stress relaxation plots for simultaneous 10% compression and twist of 10 degrees are presented in Figure 5. 
Figure 3 shows that, for the curve with the cusp, compression and shear increase linearly until they reach the level of 10% at the time of 60 seconds. Compression and shear are then maintained at the 10% level for the remainder of the 200-second simulation. It should be noted that, in the case represented by the curve with the cusp, stress reaches its maximum level at the same time that compression and shear reach their maximum levels. 
The smooth curve in Figure 3 represents the case of parabolic transition from no compression and no shear to 10% compression and 10% shear over 60 seconds, with the strain being held at that level for the remainder of the simulation. It is noted, in this case, that the maximum stress is slightly smaller (7% smaller) than it is for the linear extension regime, and that the maximum stress occurs before full compression and shear are attained. Stress relaxes by 18% of its maximum value during the period in which compression is held constant. This finding is in agreement with the standard concept that stress should not reduce to 0 in a viscoelastic solid. 
In Figure 4, four extension cases are presented. The two solid curves in this figure correspond to the two cases presented in Figure 3, in that the strain in the muscle was increased from 0 to 10% over 60 seconds in a linear (thick curve) and parabolic (thin curve) manner. The two dashed curves in Figure 4 indicate that the extension was performed over a period of 30 seconds before being maintained at the 10% level for the remainder of the 200-second simulation. As in Figure 3, the stresses depicted in Figure 4 reach their maximum levels at the same time that extension reaches its maximum level in the linear extension regime. Also, the maximum stresses are slightly smaller in the parabolic extension regime than in the linear extension regime, and these maximum stresses occur before the full extension is attained. In another finding that matches results shown in Figure 3, the stress depicted in Figure 4 relaxes by 18% of its maximum value during the period in which extension is held constant. Furthermore, Figure 4 demonstrates that as the run-up time to full extension shortens, the maximum stress increases. 
Figure 3.
Plots of normal pressure relaxation (T11) versus time for combined compression and shear of 10%, shown as a linear strain trajectory (thick line) and a parabolic strain trajectory (thin line).
Figure 3.
Plots of normal pressure relaxation (T11) versus time for combined compression and shear of 10%, shown as a linear strain trajectory (thick line) and a parabolic strain trajectory (thin line).
Figure 4.
Plots of longitudinal stress relaxation (T33) versus time for pure extension of 10%, shown as slow (thick solid line) and fast (thick dashed line) linear strain trajectories and slow (thin solid line) and fast (thin dashed line) parabolic strain trajectories.
Figure 4.
Plots of longitudinal stress relaxation (T33) versus time for pure extension of 10%, shown as slow (thick solid line) and fast (thick dashed line) linear strain trajectories and slow (thin solid line) and fast (thin dashed line) parabolic strain trajectories.
More extensive simulations were performed in the present study to represent shorter times to reach full extension in the linear regime. We hypothesized that the maximum stress value should approach the stress value previously found in the linear elastic analysis (without considering viscosity) as the time to reach full extension was decreased. That is, if extension were to take place instantaneously, we would expect there to be no time in which viscosity could act. In the simulation, we indeed found that this was the case. The maximum viscoelastic stress (T33) was 12.0 kPa when the extension was performed for the first 60 seconds of the 200-second experiment, and the stress held constant for the remainder of this experiment (Table 2). By contrast, the maximum viscoelastic stress (T33) was 18.2 kPa when the extension was performed for only the first 0.25 second of the 200-second experiment, and the stress held constant for the remainder of the experiment (Table 2). 
Table 2
Maximum Viscoelastic Stress Versus Extension Time

Extension Time, s

Maximum Viscoelastic Stress, kPa
60.012.0
30.0 13.2
5.016.5
2.0 17.0
0.25
18.2
Table 2
Maximum Viscoelastic Stress Versus Extension Time

Extension Time, s

Maximum Viscoelastic Stress, kPa
60.012.0
30.0 13.2
5.016.5
2.0 17.0
0.25
18.2
×
Figure 5.
Plots of twisting moment stress relaxation (M) versus time for simultaneous compression (10%) on the muscle element faces X1=a, X2=b and twist (10 degrees) about the muscle element axis X3, shown as a linear strain trajectory (thick line) and a parabolic strain trajectory (thin line).
Figure 5.
Plots of twisting moment stress relaxation (M) versus time for simultaneous compression (10%) on the muscle element faces X1=a, X2=b and twist (10 degrees) about the muscle element axis X3, shown as a linear strain trajectory (thick line) and a parabolic strain trajectory (thin line).
In Figure 5, twisting moment (M) relaxation is presented versus time for cases of simultaneous twist of 10 degrees and compression of 10%, each maintained at these constant values. The two curves in the figure correspond to cases in which compression and twist were increased from 0 to their maximum values over 60 seconds in a linear (thick curve) and parabolic (thin curve) manner, and then maintained at those values for the remainder of the 200-second simulations. Similar to the results in the combined compression and shear cases (Figure 3) and the pure extension cases (Figure 4), the moment in Figure 5 reaches its maximum value at the same time that compression and twist reach their maximum values in the linear extension regime. In addition, the moment relaxes by 18% of its maximum value during the time that compression and twist are held constant. Also, as in the cases presented in Figure 3 and Figure 4, the maximum moment depicted in Figure 5 is slightly smaller in the parabolic extension regime than it is in the linear extension regime, and the maximum moment occurs before the full compression and twist occur. 
Comment
The results of the present study demonstrate that because biceps muscle is viscoelastic, less force needs to be applied to it in manual therapy to produce a given deformation—provided the deformation is produced slowly (eg, 60 s)—than if the muscle were simply elastic. For quick maneuvers (eg, 0.25 s), however, this benefit of viscoelasticity is lost, and the same force is required to produce deformation as would be needed if the muscle were merely elastic. In other words, slow extension leads to lower stresses, and fast extension leads to higher stresses; so there are benefits to the patient in performing extension in a slow manner. 
Although biceps muscle tissue exhibits nonlinear behavior, the linear theory of elasticity has been used in previous studies,10 as well as in the present study, for investigating loads produced by manual therapy on biceps muscle.16 The linear theory is applicable to biceps muscle because the micro-failure region in which tissue release is felt in manual therapy occurs in the linear region.16 
The present analysis of applied loads is valid for the in vivo state and for transversely isotropic biceps muscle because the analysis uses mechanical constants that have been determined experimentally for the in vivo state, and it uses the stress-strain relation previously determined for transversely isotropic material.9,10 
In the present study, we used a parabolic trajectory for the value of the strain, E, to change from 0 to A in time T0, as indicated in equation (21). In other words, the displacement, or strain, in our model of manual therapy increases in two different trajectories of extension—linearly and parabolically. The osteopathic physician can choose to use whichever trajectory he or she prefers. 
Conclusions
In the present study, the applied loads for producing compression, shear, extension, and twist with manual therapy on biceps muscle are determined from mathematical analysis. This analysis is valid for the in vivo state and for transversely isotropic and viscoelastic biceps muscle tissue. Because the linear theory of elasticity was used in the present analysis, the results for any percentage of deformation can be obtained by linear interpolation from the results of 10% deformation provided in Table 1. 
Our results indicate that with quick maneuvers, the viscoelasticity effect is decreased (ie, 50% greater load needs to be applied to produce a given deformation). The results also indicate that 80% of the maximum stress is produced during the first 46 seconds of a 60-second OMT technique. 
The maximum stress attained using the viscoelastic model is 7% less for the same amount of deformation (ie, 10% in the present analysis) than the maximum stress attained using the elastic model. In other words, less load is needed to produce the same amount of deformation if the biceps muscle is considered as viscoelastic tissue rather than as elastic tissue. 
In a mathematical analysis of manual therapy-produced deformation of human fasciae previously published in JAOA—The Journal of the American Osteopathic Association,17 we reported that, because of the stiffness of the fasciae, no deformation was produced under physiologic loads to thick superficial fascia (eg, fascia lata and plantar fascia). However, it is logical to assume that something must yield beneath the fasciae. From the present analysis of biceps muscle, we conclude that it is the muscle beneath the fasciae that yields, and that only small loads are needed to deform this muscle. 
Because the value of Young's modulus parallel to biceps muscle fibers is approximately 15 times greater than the value perpendicular to the muscle fibers, we can conclude that the biceps muscle is 15 times stiffer in the direction parallel to the muscle fibers than in the direction perpendicular to the fibers. 
We hope that the mathematical analysis presented in the present study will strengthen the scientific foundation for the understanding of manual therapy. 
3  
Appendix Derivation of Second Terms in u1 and u2 Portions of Equation (1)
 Appendix Derivation of Second Terms in u1 and u2 Portions of Equation (1)
Appendix Derivation of Second Terms in u1 and u2 Portions of Equation (1)
 Appendix Derivation of Second Terms in u1 and u2 Portions of Equation (1) ×
We thank Ingolf Sack, PhD, of the Institute of Radiology at Universitätsmedizin Berlin, Germany, for providing information about the stiffness matrix used in this article, as well as for furnishing us with Figure 1, illustrating the plane of isotropy and the principal direction of muscle fibers in transversely isotropic biceps muscle. 
Rolf IP. Rolfing: Reestablishing the Natural Alignment and Structural Integration of the Human Body for Vitality and Well-Being. Rev ed. Rochester, Vt: Healing Arts Press;1989 .
Verala FJ, Frenk S. The organ of form: towards a theory of biological shape. J Soc Biol Struct. 1987;10:73-83.
Ward RC. Myofascial release concept. In: Basmajian J, Nyberg R, eds. Rational Manual Therapies. Baltimore, Md: Williams & Wilkins; 1993:223-241.
Plewes DB, Betty I, Urchuk SN, Soutar I. Visualizing tissue compliance with MR imaging. J Magn Reson Imaging. 1995;5:733-738.
Lewa CJ, De Certaines JD. Viscoelastic property detection by elastic displacement NMR measurements. J Magn Reson Imaging. 1996;6:652-656.
Uffmann K, Maderwald S, Ajaj W, Galban CG, Mateiescu S, Quick HH, et al. In vivo elasticity measurements of extremity skeletal muscle with MR elastography. NMR Biomed. 2004;17:181-190.
Gennisson JL, Catheline S, Chaffai S, Fink M. Transient elastography in anisotropic medium: application to the measurement of slow and fast shear wave speeds in muscles. J Acoust Soc Am. 2003;114:536-541.
Sack I, Bernarding J, Braun J. Analysis of wave patterns in MR elastography of skeletal muscle using coupled harmonic oscillator simulations. J Magn Reson Imaging. 2002;20:95-104.
Papazoglou S, Braun J, Hamhaber U, Sack I. Two-dimensional waveform analysis in MR elastography of skeletal muscles [published ahead of print March 2, 2005]. Phys Med Biol. 2005;50:1313-1325.
Lai WM, Rubin D, Krempl E. Introduction to Continuum Mechanics. Oxford, England: Butterworth-Heinemann Ltd;1996 : 303-309.
Kruse SA, Smith JA, Lawrence AJ, Dresner MA, Manduca A, Greenleaf JF, et al. Tissue characterization using magnetic resonance elastography: preliminary results. Phys Med Biol. 2000;45:1579-1590.
Ozkaya N, Nordin M. Fundamentals of Biomechanics: Equilibrium, Motion, and Deformation. 2nd ed. New York, NY: Springer; 1999: 155.
Fung YC. Biodynamics: Circulation. New York, NY: Springer-Verlag; 1984:382 .
Fung YC. A First Course in Continuum Mechanics. 3rd ed. Englewood Cliffs, NJ: Prentice Hall;1994 : 200.
Iatridis JC, Wu J, Yandow JA, Langevin HM. Subcutaneous tissue mechanical behavior is linear and viscoelastic under uniaxial tension. Connect Tissue Res. 2003;44:208-217.
Threlkeld AJ. The effects of manual therapy on connective tissues. Phys Ther. 1992;72:893-902. Available at: http://www.ptjournal.org/cgi/reprint/72/12/893. Accessed December 3, 2008.
Chaudhry H, Schleip R, Ji Z, Bukiet B, Maney M, Findley T. Three-dimensional mathematical model for deformation of human fasciae in manual therapy. J Am Osteopath Assoc. 2008;108:379-390. Available at: http://www.jaoa.org/cgi/content/full/108/8/379. Accessed December 3, 2008.
Figure 1.
Three-dimensional model of transversely isotropic biceps muscle. Transversely isotropic materials have a plane (ie, the plane of isotropy) in which the material properties are the same in different directions, but in the direction perpendicular to this plane (ie, the axis of transverse isotropy), the material properties are different than in the plane. In this model of biceps muscle, the X1, X2 plane is the plane of isotropy. The elastic properties along the X3 axis, which is the axis of transverse isotropy, are different from the elastic properties in the X1, X2 plane. (Adapted from figure provided by Ingolf Sack, PhD, of the Institute of Radiology at Universitätsmedizin Berlin, Germany, and previously published online at: http://www.elastography.de/.)
Figure 1.
Three-dimensional model of transversely isotropic biceps muscle. Transversely isotropic materials have a plane (ie, the plane of isotropy) in which the material properties are the same in different directions, but in the direction perpendicular to this plane (ie, the axis of transverse isotropy), the material properties are different than in the plane. In this model of biceps muscle, the X1, X2 plane is the plane of isotropy. The elastic properties along the X3 axis, which is the axis of transverse isotropy, are different from the elastic properties in the X1, X2 plane. (Adapted from figure provided by Ingolf Sack, PhD, of the Institute of Radiology at Universitätsmedizin Berlin, Germany, and previously published online at: http://www.elastography.de/.)
Figure 2.
Three-dimensional model of skeletal muscle element used in analysis of displacement caused by applied loads. OX1, OX2, and OX3 are the axes in the Cartesian coordinates system; F denotes the longitudinal stress; M denotes the moment; N and T denote the normal pressure and tangential force, respectively. The six faces of this model are designated by combinations of the capital letters A, B, C, D, E, G, H, and I. The symbols a, b, -b, and c represent the coordinate values of the faces along the three axes.
Figure 2.
Three-dimensional model of skeletal muscle element used in analysis of displacement caused by applied loads. OX1, OX2, and OX3 are the axes in the Cartesian coordinates system; F denotes the longitudinal stress; M denotes the moment; N and T denote the normal pressure and tangential force, respectively. The six faces of this model are designated by combinations of the capital letters A, B, C, D, E, G, H, and I. The symbols a, b, -b, and c represent the coordinate values of the faces along the three axes.
Figure 3.
Plots of normal pressure relaxation (T11) versus time for combined compression and shear of 10%, shown as a linear strain trajectory (thick line) and a parabolic strain trajectory (thin line).
Figure 3.
Plots of normal pressure relaxation (T11) versus time for combined compression and shear of 10%, shown as a linear strain trajectory (thick line) and a parabolic strain trajectory (thin line).
Figure 4.
Plots of longitudinal stress relaxation (T33) versus time for pure extension of 10%, shown as slow (thick solid line) and fast (thick dashed line) linear strain trajectories and slow (thin solid line) and fast (thin dashed line) parabolic strain trajectories.
Figure 4.
Plots of longitudinal stress relaxation (T33) versus time for pure extension of 10%, shown as slow (thick solid line) and fast (thick dashed line) linear strain trajectories and slow (thin solid line) and fast (thin dashed line) parabolic strain trajectories.
Figure 5.
Plots of twisting moment stress relaxation (M) versus time for simultaneous compression (10%) on the muscle element faces X1=a, X2=b and twist (10 degrees) about the muscle element axis X3, shown as a linear strain trajectory (thick line) and a parabolic strain trajectory (thin line).
Figure 5.
Plots of twisting moment stress relaxation (M) versus time for simultaneous compression (10%) on the muscle element faces X1=a, X2=b and twist (10 degrees) about the muscle element axis X3, shown as a linear strain trajectory (thick line) and a parabolic strain trajectory (thin line).
Table 1
Stress Values for Applied Forces in Elastic and Viscoelastic Models


Stress Value
Applied Force*
Elastic Model
Viscoelastic Model
Compression and Shear, kPa (lb/in2) (k1=0.1; k3=0.9; k5=0)
□ Compressive stress -1.30 (-0.19) -0.80 (-0.12)
□ Shear stress10.80 (1.57)6.64 (0.96)
Longitudinal Extension, kPa (lb/in2) (k1=0; k3=1.0; k4=0.1; k5=0)
□ Longitudinal stress18.50 (2.68)11.40 (1.65)
Twist and Compression, kPa-cm§ (k1=0; k2=0.9; k3=0.9; ψ=10 degrees)
□ Twisting moment23.7014.60
 Abbreviations: k1 denotes shear ratio resulting from application of tangential force; k2, extension ratio resulting from compression and shear on surface of muscle; k3, compression ratio resulting from applied normal pressure; k4, extension ratio resulting from applied longitudinal force; k5, Cosψ - 1; ψ, angle caused by twist applied to muscle.
 *Applied 10% compression and applied 10% shear are on face X1=a of skeletal muscle element. Applied extension is parallel to fiber direction on face X3=c of skeletal muscle element. See Figure 2.
 k2 and k4 can be determined from equation (6) in text.
 k2 can be determined from equation (9) in text.
 §k4 can be determined from equation (11) in text.
 Twisting moment per unit area about the X3 axis.
Table 1
Stress Values for Applied Forces in Elastic and Viscoelastic Models


Stress Value
Applied Force*
Elastic Model
Viscoelastic Model
Compression and Shear, kPa (lb/in2) (k1=0.1; k3=0.9; k5=0)
□ Compressive stress -1.30 (-0.19) -0.80 (-0.12)
□ Shear stress10.80 (1.57)6.64 (0.96)
Longitudinal Extension, kPa (lb/in2) (k1=0; k3=1.0; k4=0.1; k5=0)
□ Longitudinal stress18.50 (2.68)11.40 (1.65)
Twist and Compression, kPa-cm§ (k1=0; k2=0.9; k3=0.9; ψ=10 degrees)
□ Twisting moment23.7014.60
 Abbreviations: k1 denotes shear ratio resulting from application of tangential force; k2, extension ratio resulting from compression and shear on surface of muscle; k3, compression ratio resulting from applied normal pressure; k4, extension ratio resulting from applied longitudinal force; k5, Cosψ - 1; ψ, angle caused by twist applied to muscle.
 *Applied 10% compression and applied 10% shear are on face X1=a of skeletal muscle element. Applied extension is parallel to fiber direction on face X3=c of skeletal muscle element. See Figure 2.
 k2 and k4 can be determined from equation (6) in text.
 k2 can be determined from equation (9) in text.
 §k4 can be determined from equation (11) in text.
 Twisting moment per unit area about the X3 axis.
×
Table 2
Maximum Viscoelastic Stress Versus Extension Time

Extension Time, s

Maximum Viscoelastic Stress, kPa
60.012.0
30.0 13.2
5.016.5
2.0 17.0
0.25
18.2
Table 2
Maximum Viscoelastic Stress Versus Extension Time

Extension Time, s

Maximum Viscoelastic Stress, kPa
60.012.0
30.0 13.2
5.016.5
2.0 17.0
0.25
18.2
×
Appendix Derivation of Second Terms in u1 and u2 Portions of Equation (1)
 Appendix Derivation of Second Terms in u1 and u2 Portions of Equation (1)
Appendix Derivation of Second Terms in u1 and u2 Portions of Equation (1)
 Appendix Derivation of Second Terms in u1 and u2 Portions of Equation (1) ×