Dynamic Model of the Octopus Arm. I. Biomechanics of the Octopus Reaching Movement
Abstract
The octopus arm requires special motor control schemes because it consists almost entirely of muscles and lacks a rigid skeletal support. Here we present a 2D dynamic model of the octopus arm to explore possible strategies of movement control in this muscular hydrostat. The arm is modeled as a multisegment structure, each segment containing longitudinal and transverse muscles and maintaining a constant volume, a prominent feature of muscular hydrostats. The input to the model is the degree of activation of each of its muscles. The model includes the external forces of gravity, buoyancy, and water drag forces (experimentally estimated here). It also includes the internal forces generated by the arm muscles and the forces responsible for maintaining a constant volume. Using this dynamic model to investigate the octopus reaching movement and to explore the mechanisms of bend propagation that characterize this movement, we found the following. 1) A simple command producing a wave of muscle activation moving at a constant velocity is sufficient to replicate the natural reaching movements with similar kinematic features. 2) The biomechanical mechanism that produces the reaching movement is a stiffening wave of muscle contraction that pushes a bend forward along the arm. 3) The perpendicular drag coefficient for an octopus arm is nearly 50 times larger than the tangential drag coefficient. During a reaching movement, only a small portion of the arm is oriented perpendicular to the direction of movement, thus minimizing the drag force.
INTRODUCTION
The biomechanical principles of movement generation differ in animals with and without a rigid skeleton. The latter possess a hydrostatic skeleton, consisting mainly of muscles and fluid, which because they are incompressible, support the flexible body. Hydrostatic skeletons are generally of 2 types. In muscular hydrostats, muscles and other tissues form a solid structure without a separate enclosed fluid volume (e.g., cephalopod tentacles, elephant trunk, vertebrate tongue; Kier and Smith 1985). In a second type of hydrostatic skeleton, muscles compose a body wall that surrounds a fluidfilled cavity (e.g., sea anemones and worms, Kier 1992).
All muscular hydrostats consist of closely packed arrays of muscle fibers organized in 3 main directions (Kier and Smith 1985): parallel, perpendicular, and helical or oblique to the long axis. The major feature of muscular hydrostats is that they maintain a constant volume (Kier and Smith 1985). This constraint allows forces to be transferred between the longitudinal and the transverse directions. All motion in a muscular hydrostat is based on combinations of 4 elementary movements that can occur at any location: elongation, shortening, torsion, and bending (Kier and Smith 1985). Muscular hydrostats thus have an exciting and unique property: the muscles not only generate the forces required to implement movements, but also provide the skeletal support, i.e., the muscular system can function as a dynamic skeleton. This allows smooth changes of shape, resulting in a potentially vast repertoire of movements.
While an articulated arm has a limited number of degrees of freedom (DOFs), the octopus arm has virtually an infinite number of degrees of freedom and thus poses a considerable challenge from the control perspective.
The uniqueness of this system raises fundamental questions: What constraints does the mechanical structure of the muscular hydrostat structure impose on the motor control system? How does the motor control system select a specific movement from an apparently infinite number of available possibilities? What does the movement repertoire of the octopus arm teach us about evolutionary and optimization processes and the constraints that have shaped the octopus motor control system? How does the nervous system control and coordinate the muscles?
Understanding the principles that govern the motion of this arm is not only interesting from the biological standpoint but is also important for robotics engineering. Researchers in this field face similarly complex problems when designing control systems for flexible robotic arms with a large number of DOFs (e.g., Chirikjian and Burdick 1994; Takanashi et al. 1993; Walker 2000).
Reaching movement
The octopus uses a limited number of stereotypical motion patterns to reach out to catch and inspect an object (Gutfreund et al. 1996, 1998). A bend is created, usually near the base of the arm (but it can be created almost anywhere along the arm), and is propagated along the arm toward the tip. The arm segment proximal to the bend is maintained relatively straight and the bend is always curved dorsally (Fig. 6A). The path of this bend is confined to a plane and the bend propagates according to an invariant basic velocity profile. Gutfreund et al. (1996) proposed that this mode of reaching solves the problem of multiple DOFs by reducing their potentially large number to a very small number of control parameters that shape the movement—one DOF for the motion of the bend along the arm and 2 additional DOFs for the orientation of the base of the arm in space. Sumbre et al. (2001) suggested that the brain issues commands to orient the arm direction and set the extension parameters, whereas the details for muscle activation are represented and carried out within the neuromuscular system of the arm itself.
Bending results from a selective contraction of longitudinal muscles on one side of the arm together with the contraction of transverse muscles to resist shortening (Kier 1992; Kier and Smith 1985). The pattern of muscle activation that created the bend can then be propagated along the arm creating an extension movement. Gutfreund et al. (1998) proposed a different mechanism for bend propagation—a stiffening wave caused by a symmetrical muscle activation pattern that propagates along the arm and propels an existing bend. The model of the octopus arm that is presented here was developed to test the plausibility of the stiffening wave hypothesis and how such a wave can be controlled to produce extension movements.
Modeling muscular hydrostats
We do not fully understand the constraints imposed by the biomechanics of muscular hydrostats on their control systems, nor the strategies these control systems adopt. Because these cannot be studied only by an experimental approach, several mathematical models have been developed to characterize muscular hydrostats. Most of these muscular hydrostat models are not sufficiently general to describe the motion of the octopus arm for one or more of the following reasons.
The models are static or quasistatic, i.e., they do not describe the full dynamics of motion: Wadepuhl and Beyn (1989) developed a quasistatic model of wormlike forms; Skierczynski et al. (1996) and Kristan et al. (2000) developed a quasistatic model of the hydrostatic skeleton of the medicinal leech; and Wilson et al. (1991) gave a static continuum model of the elephant trunk. Related to these models are kinematic models; Drushel et al. (1998) constructed 2 kinematic models of the buccal mass of Aplysia californica. (See Drushel et al. 2002 and Neustadter et al. 2002 for the development of these models.)
The models are restricted to describe specific movements characterized by a low number of DOFs. Van Leeuwen and Kier (1997) developed a dynamical model of the squid tentacles capable of producing fast elongation extension movements. This model can describe only the forward extension of the tentacle.
The models cannot account for all feasible motions (e.g., they assume some constraints on the shape of the structure). Chiel et al. (1992) produced a quantitative model of the reptilian tongue, but they neglected mass and inertial forces, which are likely to be small compared to the muscle forces involved in the movements they modeled.
The models do not consider external forces like gravity or water drag forces, which need to be considered when modeling unconstrained movements of the octopus arm (e.g., Alscher and Beyn 1998; Wadepuhl and Beyn 1989).
The models do not incorporate the constant volume constraint. For example, Jordan (1996) developed a dynamical model of the leech that couples internal factors (masses, muscle mechanics, soft tissue mechanics, and muscle activation patterns) and external factors (water hydrodynamics) to predict the swimming behavior. This model used a rigid endoskeleton as a first approximation to the hydrostatic skeleton.
Some of these problems are remedied by our dynamic model of the octopus arm, based on physiological and kinematic studies. This model is currently essentially 2dimensional (2D), although it allows natural extension to 3 dimensions.
Here we show that the physics of the arm and its interaction with water play a significant role in shaping the arm's movements. We show that an extension movement similar to the natural reaching movement can be produced by a simple neural command propagating along the arm. This produces a wave of muscle contraction that travels along the arm, causing it to extend. In paper II, the second part of this zpart series (Yekutieli et al. 2005), we investigate the control of reaching movements, focusing on the properties of the command that creates the wave of muscle contraction. We show a relationship between the amplitude and timing of this command, which may be used by the motor control system to choose the minimal forces needed to extend an arm according to a desired kinematic plan. Our model predicts the possibility of using higher forces to resist external perturbations while keeping the same kinematics. This prediction suggests that the octopus may be able to independently control 2 major aspects of arm extension: the kinematics of the movement and its resistance to perturbations.
A preliminary report on our work was presented in Aharonov et al. (1997).
THE MECHANICAL MODEL—AN OVERVIEW
The arm is modeled as a 2D array of point masses and springs. All the arm's mass is concentrated in discrete point masses connected by massless idealized muscles that are modeled as damped (linear or nonlinear) springs (Fig. 1). The 2n masses are arranged in n pairs, each consisting of one ventral and one dorsal mass. The longitudinal muscles and the transverse muscles connect the masses. Through this arrangement, every 2 adjacent pairs of masses and their corresponding muscles enclose a quadrilateral compartment. These (n − 1) compartments play a crucial role in the model because the constant volume constraint is enforced locally in each of the compartments. [See Fig. 10 in Kier (1988) for the arrangement of the musculature of the octopus arm.]
Our model is 2D in the sense that all force vectors lie within the x–y plane, where −y is the direction of gravity; thus the arm's motion is constrained to be planar. The following notations are used throughout:
 x
 a scalar
 r⃗
 a vector in ℜ^{2} or ℜ^{3} of a physical quantity like velocity
 r
 a vector in ℜ^{m}, m > 3
 x̄
 a scalar average value
 X
 a matrix
 x_{i}, x_{ij}
 elements of a vector and matrix, respectively
 ẋ, ẍ
 time derivatives of x and ẋ, respectively
Model assumptions and limitations

The movements that the model can simulate are free movements that do not interact with objects or the seabed.

A discrete, rather than a continuous description of the arm is used to simplify modeling. There are 20 segments in the model, and this is assumed a sufficiently accurate approximation to the continuous structure of the octopus arm; modeling an arm with 40 segments made little difference (see results).

All forces included in the model are confined to one plane, which limits its generality but significantly reduces the computational cost. Moreover, the reaching movement that is studied here is mainly planar.
THE MODEL IN DETAIL
Arm properties
The physical dimensions of the model arm approximate those of a typical Octopus vulgaris arm in its relaxed state, i.e. when it is submerged in seawater and not subjected to external forces apart from its own weight. Arm width is 15 mm at the base and 1 mm at the tip. Arm length is 240 mm and the depth at every point equals the width at that point and remains constant throughout the simulation (Fig. 1).
The model forces
Four types of forces act on the arm's masses:
The internal forces generated by the arm muscles (f^{m}).
The vertical forces resulting from the combined influences of gravity and arm buoyancy in seawater (Archimedes' law of buoyancy) (f^{g}).
The drag produced by arm motion through the water (f^{w}).
The internal forces responsible for maintaining the constant volume constraints (f^{c}).
Model input
The only input to the model is an activation pattern representing the neural input to the muscular system (described below). The initial shape of the arm is a typical shape that real octopus arms assume before a reaching movement.
Equations of motion and numerical solution
The motion equations (Newton's second law) can be compactly written as (1) where M is a diagonal mass matrix and q is the position vector composed of the coordinates of the 2n masses (see appendix A for the definitions of matrices and vectors). This equation is numerically integrated using the explicit Runge–Kutta (2,3) method with an adaptive step size (Bogacki and Shampine 1989). The initial conditions are the initial positions of all the masses. Initial velocities are set to zero. The simulations were written in Matlab version 6.5 (The MathWorks).
The arm's weight
According to Archimedes' law of buoyancy, the arm's weight in water can be expressed as (2) where ρ is density, 1,022 kg/m^{3} for seawater and 1,042 kg/m^{3} for an octopus arm; V_{arm} is a 4n × 4n matrix of segment volumes; and g is the gravitational acceleration of the Earth.
The total weight f^{g} in Eq. 2 is distributed among the arm's compartments according to their volume. Within each compartment, the weight is divided equally among the 4 masses, taking into account that all but the first and last mass pairs are shared by 2 neighboring compartments.
The constraint forces
The constant volume constraint for each compartment is maintained by pressure changes within the compartment. This pressure induces forces exactly counteracting those that would otherwise cause volume changes. The transformation from the compartmental pressures p into the constraint forces f^{c} is linear and can be written compactly as (3) Note that, by using the form above, we reduce the dimension of the problem from 2n (number of point masses) to n − 1 (number of the compartmental pressures). The matrix C is found using the Lagrange formulation with Lagrange multipliers corresponding to the compartmental pressures, similar to Wadepuhl and Beyn (1989) (see appendix B).
Equation 1 cannot be solved explicitly because f^{c} cannot be directly determined. The constraint force f^{c} must be derived indirectly from the constant volume constraints (see appendix B). The result of the derivation is the pressure vector p (4) (See appendix B for the definition and derivation of the matrix G and the vector γ.)
Once the compartmental pressures p are known, they can be linearly transformed into f^{c} (Eq. 3). All forces are now expressed explicitly and the equations of motion can be solved (5)
The drag forces
A body moving through fluid is subjected to drag forces that increase with the relative velocity of the body. This dependency of drag force on velocity is rather complex because it is composed of a direct dependency on the squared velocity, and an indirect effect through the Reynolds number,^{1} which itself depends on velocity. The characteristic velocities of octopus arms during reaching movements in seawater lie in the range of 20–60 cm/s (Gutfreund et al. 1996). The Reynolds number of an octopus arm for these conditions is ≅3000 (see details in calculations and simplifications below). A value of Re between 50 and 200,000 indicates a moderately turbulent flow.
When the arm moves in water, each of the segments moves with some velocity v that can be decomposed into 2 components: v_{per} in the direction perpendicular to the long axis of the segment and v_{tan} in the direction tangential to the long axis. Similarly, the water drag forces acting on the segment can be decomposed into a perpendicular component d_{per} and a tangential component d_{tan}.
Both the perpendicular and the tangential components of the drag forces are calculated for each segment and then decomposed into their x and y components. Following Newton's second law and the dimensional analysis technique,^{2} the two drag components in steady flow^{3}are (Vogel 1981) (6) (7) where P_{a} is the projected area of the segment in the plane perpendicular to the direction of ν⃗_{per} and S_{a} is the surface area of the segment. The components of the segment's velocity in the perpendicular and tangential directions are ν⃗_{per} and ν⃗_{tan}, respectively. The coefficients c_{per} and c_{tan} express the perpendicular and tangential drag coefficients, respectively, which are functions of Reynolds number as follows The functions f and g are of the form u_{1} + u_{2}(Re)^{u3}, where u_{1}, u_{2}, and u_{3} are scalars (Vogel 1981).
A biologically inspired model for swimming in the chaetognath Sagitta elegans (Jordan 1992) was based on the drag force equations (Eq. 6 and 7). It used a perpendicular drag coefficient (c_{per}) from curves fitted for cylinders in cross flow, ^{4} and a tangential drag coefficient (c_{tan}) from curves fitted for flat plates in streamwise flow.^{5} The functions expressing the perpendicular and tangential drag coefficients were taken from Vogel (1981) and White (1974) and were expressed as follows where Re_{per}, the Reynolds number for the perpendicular direction, is (2rv)/υ, where r is the cylinder radius and υ is the kinematic viscosity. Reynolds number for the tangential direction Re_{tan} is (xv)/υ, where x is the cumulative cylinder's length.
Experimental estimation of the drag forces
We did not use the above tangential and perpendicular drag coefficients (Jordan 1992) but instead experimentally estimated the drag forces for the following reasons: 1) The drag coefficients depend strongly on the structure and surface properties of the moving body. The shape of the octopus arm is very complex, but Jordan's values were obtained for idealized shapes. 2) It is essential to use as accurate values as possible because we found that the behavior produced by the model is sensitive to changes in the values of the drag coefficients (especially the tangential one; see below in results).
An experiment was used to evaluate the 2 components of the drag force, d_{per} and d_{tan}. First we calculated the real projected area P_{a} and sectional area S_{a}. We then used Eq. 6 and 7 to obtain the experimental drag coefficients c_{per} and c_{tan}. These values were inserted into the model as described next.
We used a longitudinal segment of an amputated octopus' arm (for details on surgical procedure, see Sumbre et al. 2001). A thin metal rod (weight 0.2 g) was inserted into the segment to maintain its shape by eliminating the natural flexibility of the arm. The dimensions of the segment were: weight (including the metal rod) 5.29 g, volume 5 ml, and length 82 mm. The radius of the distal tip of the segment was 2.5 mm and the radius of its proximal tip was 7 mm.
The segment was attached to different metal loads by a nylon string of negligible mass and was dropped through seawater in a 75 × 35 × 40cm (length × width × depth) aquarium. We dropped the segment so as to let it fall horizontally and vertically. For the former (Fig. 2A), a horizontally oriented arm segment was pulled down by 8 different loads between 1.688 to 7.26 g. Here the velocity of the segment is confined to the direction perpendicular to the long axis of the segment. In the vertical fall (Fig. 2B), a vertically oriented arm segment was pulled down by 13 different loads between 0.145 to 0.947 g. The segment's velocity was confined to the direction tangential to the long axis of the segment. All falls were repeated 4–6 times for each load. Images of the falling arm segment were recorded by a PAL SVHS video camera with 40 ms between adjacent frames.
CALCULATIONS AND SIMPLIFICATIONS.
When the segment is dropped, it accelerates until its velocity levels off at some constant value, i.e., no net force acts on the segment (by Newton's second law where f⃗ are force vectors). In our experiment, the relevant forces acting on the segment were confined to the vertical direction.
The force equations are (8) (9) where f_{Ase} = V_{se}ρg⃗ is Archimedes' force acting on the segment and f_{Alo} = V_{lo}ρg⃗ is Archimedes' force acting on the load. The volumes of the segment and the load are V_{se} and V_{lo}, respectively, and the weight of the segment and the load are W_{se} and W_{lo}.
In Eq. 8 and 9 all forces except the drag forces are known or can be easily determined, so the extraction of the drag forces is trivial.
To estimate the projected area p and the surface area s, we used a truncated cone as a simplification for the more complex structure of the octopus arm. The projected area P_{a} is 2r̄l, where r̄ is the mean segment radius of the arm and l is the segment length. The sectional surface area S_{a} is 2πr̄l.
For the perpendicular component Re = 2r̄v̄_{per}/υ. For the tangential component Re = lv̄_{tan}/υ, where l is the segment length, v̄ is the segment average velocity, r̄ is the segment average radius, and υ is the kinematic viscosity (υ = 1 × 10^{−6} m^{2}/s).
ANALYSIS OF THE EXPERIMENT.
For each experiment, we plotted velocity versus time, found the levelingoff velocity, and computed its average for the particular load. We then plotted the average velocities versus the different loads. Figure 3A shows the perpendicular drag components versus squared velocity as obtained from the horizontal fall experiments; Fig. 3B gives the tangential drag components from the vertical fall experiments. A linear fit was sufficient and we used the slope to calculate the drag coefficients.
For the perpendicular component, the slope of the linear fit is 0.5ρP_{a}c_{per} (Fig. 3A and Eq. 6), and for the tangential component the slope of the linear fit is 0.5ρS_{a}c_{tan} (Fig. 3B and Eq. 7).
EXPERIMENTAL RESULTS.
Drag coefficients are dependent on Reynolds number but, because of the linear relations seen in Fig. 3, we can simplify the calculations by regarding the drag coefficients as constants. For the perpendicular drag coefficient c_{per} (from Eq. 6)
For the tangential drag coefficient c_{tan} (from Eq. 7)
Note that the perpendicular drag coefficient is nearly 50fold larger than the tangential drag coefficient. This difference may be one of the factors that shaped the octopus reaching movement such that only a small portion of the arm is oriented perpendicular to the direction of movement, thus minimizing the drag force (see discussion).
Muscle forces
We used 2 types of muscle models.
A nonlinear model similar to the models described by Van Leeuwen and Kier (1997) and Zajac (1989). It includes nonlinear force–length and force–velocity relations.
A linear damped spring model (similar to the model of Ekeberg 1993).
Most of our work used the nonlinear muscle model. This is compared with the linear model below.
We assume that the arm's muscle fibers are homogeneous and evenly distributed throughout the arm and that all muscles have similar functional properties (Matzner et al. 2000). The various parameters of each muscle model are scaled for the dimensions of the different muscles according to the dimensions of a real octopus arm in a relaxed state.
For simplicity, we present only the magnitude of the force and not its direction, but all real calculations are vectorial, taking the direction into account.
The nonlinear muscle model
Every muscle in the model exerts the following force (10) where A is the crosssectional area of the muscle, a ∈ [0, 1] is a dimensionless activation function, l(t) is muscle length, and l_{0}^{m} is the specific length of the muscle at which active muscle force achieves its peak value (terminology after Zajac 1989). The active stress F_{a} equals the normalized active force (Fig. 4A) multiplied by the maximum stress. The dimensionless normalized force–velocity relation is expressed by F_{v}, which is a function of v/v_{max}, where v = dl(t)/dt, the rate of change of muscle length with time, and v_{max} is the absolute value of the maximum shortening velocity. The passive stress F_{p} is equal to the normalized passive force (Fig. 4A) multiplied by the maximum stress.
The output force is controlled by the activation function a(t) whose range is [0, 1]. When the activation is zero, the muscle is passive and acts as a nonlinear spring according to F_{p} (Fig. 4A, “passive” curve). In a fully activated isometric contraction the muscle force is the sum of the active and passive contributions F_{a} and F_{p} (Fig. 4A, “total” curve). When muscle length changes, its output force is also influenced by the rate of change of muscle length. This nonlinear relation is introduced by F_{v} (see Fig. 4B and explanation below).
CHOICE OF PARAMETERS IN THE NONLINEAR MUSCLE MODEL.
The relaxed crosssectional area differs between the longitudinal and transverse muscles (Fig. 1) where w_{r} and l_{r} are the relaxed lengths^{6} of the transverse and longitudinal segmental muscles, respectively; and d is the depth of the segment, which equals w_{r} at the initial relaxed shape of the arm. The factor of in A^{longitudinal} is ascribed to the fact that there are a total of 2n longitudinal muscles in the model for n segments (and only n + 1 ≅ n transverse muscles in the model).
The activation dependency on time is expressed by a(t) and described below. Note that we assume a linear scaling of the force with the activation parameter. This is an approximation to experimental data relating stimulus strength and frequency to muscle force in muscles of other cephalopods (squid and sepia mantle muscles; Milligan et al. 1997, Figs. 2 and 7).
Values of F_{a}, F_{p}, and l_{0}^{m}: The active and passive muscle force curves (Fig. 4A) are based on spline smoothing of experimental data from a lateral group of longitudinal muscles from Octopus vulgaris arm. The curves were obtained from isometric measurements of muscle force with and without muscle activation (H. Matzner and B. Hochner, unpublished results). The general shape of the curves is similar to previously reported active and passive force–length relations (e.g., funnel retractor muscle of Octopus vulgaris; Lowy and Millman, 1962, Fig. 4; obliquely striated muscle from cuttlefish mantle; Milligan et al. 1997, Fig. 5B; smooth muscle I2 in the buccal mass of Aplysia; Yu et al. 1999, Fig. 2C).
The peak stress in a brief isometric tetanus was found to be 2.62 × 10^{5} N/m^{2} in squid mantle muscle and 2.26 × 10^{5} N/m^{2} in cuttlefish (Curtin et al. 2000; Milligan et al. 1997). Peak stress in a brief isometric tetanus was 1.31 × 10^{5} N/m^{2} in squid tentacle muscles and 4.68 × 10^{5} N/m^{2} in the squid arm (Kier and Curtin 2002). The maximum value observed from twitches of the funnel retractor muscles of octopus was 5.0 × 10^{5} N/m^{2} (Lowy 1962). Our measurements on a lateral group of longitudinal octopus arm muscles gave 4.0 × 10^{4} N/m^{2} (H. Matzner and B. Hochner, unpublished results). The maximal active muscle stress used in the model here is 8.0 × 10^{4} N/m^{2}.
Both the active and passive force–length curves are normalized to l_{0}^{m}, the length at which active muscle force is maximal. The relaxed length of the group of muscles we measured was 0.8l_{0}^{m} (H. Matzner and B. Hochner, unpublished results), similar to that reported for octopus funnel retractor muscle of Lowy and Millman (1962), Fig. 4) for Sepia mantle muscles (Milligan et al. 1997, Fig. 5A) and for smooth muscle I2 in the buccal mass of Aplysia (Yu et al. 1999, Fig. 2C).
The force–velocity properties (Fig. 4B), which are not known for the octopus arm muscles, were modeled using apiecewise linear approximation to the force–velocity relations described by Zajac (1989). In our approximation, the numerical values of F_{v} at the boundaries are the same as those in the relation described by Zajac (1989). We use the following relation where v* = v/v_{max} is the normalized velocity and v_{max} is the maximum shortening velocity taken as 10 muscle lengths s^{−1} (Zajac 1989). This value has been found in many systems and our preliminary measurements comply with it. Kier and Curtin (2002) reported a maximum shortening velocity of 15.5l_{0}^{m} s^{−1} for squid tentacle fibers and 1.5l_{0}^{m} s^{−1} for squid arm muscle fibers.
The linear muscle model
Every linear muscle in the model exerts the following force (11) where l_{rest} is the rest length of the muscle. This was chosen as the largest length at which both active and passive forces are zero in real muscles. The linear damping coefficient α has dimensions of Ns/m. The passive spring constant of the muscle is expressed by k^{0}, and the maximal active spring constant of the muscle by k^{max}, both have dimensions of N/m.
PARAMETER CHOICE FOR THE LINEAR MUSCLE MODEL.
Values of k_{max}, k_{0}, and l_{rest}: The various parameters of the linear muscle model were chosen by approximating their counterparts in themore realistic nonlinear model as depicted in Fig. 4C. The spring rest length l_{rest} is the length at which both models produce no force. Its value is 0.4l_{0}^{m} (0.4 of the length at which the active muscle force peaks in the nonlinear muscle model).
Hooke's law can also be expressed by the ratio of stress to strain (length change divided by original length). This ratio is termed stiffness (Curtin et al. 2000) and is used here in the calculation of the spring constants. The total stiffness T_{total} is the sum of the passive stiffness T_{passive} and the active stiffness T_{active}. It was taken as the slope of the line connecting the point (l_{rest}, 0) and the maximal force point (1, 1) in Fig. 4C. This gave a value of T_{total} as 1.34 × 10^{5} N/m^{2}.
The approximation of the measured passive force–length relationship (Fig. 4C) by a linear relationship is not trivial because the measured relationship is highly nonlinear. Using the nonlinear model as a reference, we searched for a value for the linear passive force that best fits this reference. We ran 2 simulations without any active force, using either the linear or nonlinear muscle models, with different linear passive values for every run, until the shape changes of the arm in the linear case matched those of the nonlinear case. The optimal value for T_{passive} was found to be 2,000 N/m^{2}. The active stiffness T_{active} is the difference between T_{total} and T_{passive} and its value was 1.32 × 10^{5} N/m^{2}. The maximal spring constant k^{max} equals T_{active} times the muscle's relaxed crosssectional area, divided by its relaxed length. Thus, in Eq. 11 the calculation of k^{max} is different for a longitudinal and a transverse muscle.
For longitudinal muscles whereas for the transverse muscles The passive component of the muscle force (k^{0} in Eq. 11) is similarly computed using the value of T_{passive} instead of T_{active}. The property α of Eq. 11 is calculated analogously to k^{max}, where T_{active} is substituted by α_{0} = 9 Ns/m^{2}, which is also assumed to be a constant property of all muscle fibers. This value was chosen such that the damping forces produced by the linear model are of the same order of magnitude as the forces produced by the force–velocity part of the nonlinear muscle model.
CHOOSING THE ACTIVATION FUNCTION FOR BOTH MUSCLE MODELS.
The general activation scheme used here is a uniform stiffening wave that equally activates all muscles in a single segment and propagates along the arm from base to tip. Constraining all the muscles in each segment to be equally activated reduces the number of control parameters from 3n − 2 to n − 1, thus significantly reducing the complexity of the control problem. As we show below, this reduction in complexity does not hinder the generation of the desired arm movements. We used either a constant velocity or a bellshaped velocity profile for the propagation of the activation signal.
In general, for a constant velocity profile the activation of the muscles can be expressed as (12) where ac_{a} is the maximum value of the activation parameter and lies between 0 and 1; i is the segment index of the muscle i = 1 . . . n − 1, where i = 1 is the most proximal segment; τ (s/segment) is the time for the signal f to pass one segment; and i_{0} is the phaseshift parameter that is equal for all segments and was set at i_{0} = 1. Once the activation function f is chosen, the only remaining control parameters are ac_{a}, i_{0}, and τ.
One such activation function is an inverted sigmoid (Fig. 5) (13) This sigmoid shape of the activationtime curve (Fig. 5A) approximates the overall shape of the EMG activation measured in an octopus arm during reaching movements (Gutfreund et al. 1998). The sustained muscle activity (amplitude of the EMG signal distal to the bend in the arm) was lower than the maximal level. We tried a pattern of activation with lower sustained activity in the model and it produced results similar to those using the sigmoid shape.
The time taken by the activation wave front to pass along the whole arm or activation traveling time is marked by ac_{t}. It equals τ (the time the signal takes to pass one segment) multiplied by the number of segments. Note that in this model there is no delay between the activation of a muscle and the consequent contraction. Assuming the excitation–contraction delay is constant (Gutfreund et al. 1998) and that there is a onetoone relation between excitation and contraction, this delay may be absorbed into the phaseshift parameter i_{0}.
We used a minimum jerk polynomial^{7} for the bellshaped form of the signal velocity. The maximal signal velocity was scaled such that the signal traveled along the arm from base to tip in a given time with the desired velocity (i.e., the integral of the velocity over time equals the length of the arm). Thus for a given total movement time, the average signal velocity was the same as the constant velocity (Fig 5C).
RESULTS
Drag force
During a reaching movement, the proximal segment of the arm (the part between the bend and the base of the arm) does not move and therefore does not suffer from drag forces. The area around the bend is relatively small but is exposed to large drag forces reflecting the large perpendicular coefficient. In contrast, the distal part (from the bend to the tip of the arm) has on average a much larger surface area, although its drag is smaller because of the nearly 50 times smaller tangential coefficient. Measuring the total tangential and perpendicular drag forces integrated over the entire simulated reaching movement revealed that they have almost the same magnitude.
Simulations and experimental results
We now describe behaviors observed in octopuses and the results of the simulations mimicking such behaviors. In the simulations, we change the input to the model (analogous to a neural command) by changing the activation amplitude ac_{a} and the activation traveling time ac_{t}, while keeping all other parameters fixed. Unless otherwise stated, we use the nonlinear muscle model.
Producing a reaching motion
OBSERVED BEHAVIOR.
It is easy to induce an octopus to make stereotypical reaching movements by holding a target in front of it (Gutfreund et al. 1996). The octopus extends one or more of its arms toward the target using a wavelike propagation of a bend toward the tip, while keeping the proximal part of the arm almost straight (Fig. 6A). The bend travels along a simple planar, slightly curved path, connecting the animal's body and the target. Normalized tracings of the position of the point of maximum curvature as it progresses along the arm in movements of varying durations and magnitudes reveal a typical invariant shape for the point's tangential velocity profile (Gutfreund et al. 1996). A denervated arm produces a similar movement when electrically stimulated (Fig. 6B; Sumbre et al. 2001).
SIMULATION RESULTS.
The model is capable of reproducing the empirical characteristics of natural bend propagation. Moreover, this seemingly complex motion can be reproduced by a very simple muscle activation plan, consisting of a uniform stiffening wave propagating along the arm from the base to the tip at a constant velocity. Figures 6, C and D and 7A show traces from simulated extension movements. Replicating a typical position assumed by the octopus at the onset of reaching, the arm is curled dorsally near the animal's body with its ventral surface up. Then the activation wave front progresses distally activating more and more segments. As in natural and evoked extension movements, the part of the arm from the base to the bend is almost straight.
Bend point path and velocity profile
OBSERVED BEHAVIOR.
Octopus reaching movements are characterized by an invariant tangential velocity^{8} profile of the bend point (Gutfreund et al. 1996). This profile can be divided into 3 phases: 1) a monotonic decrease in velocity or by an initial increase in velocity followed by a local minimum; 2) a monotonic increase in velocity; and 3) a decrease in the velocity until the bend point disappears (Fig. 8A). Velocity profiles of denervated arms show similar patterns (Fig. 8B).
SIMULATION RESULTS.
The shape of the simulated arm was displayed on a screen and the position of the bend point was marked manually, using the same method as in Gutfreund et al. (1996). The path of the bend point (Fig. 7B) was an almost straight line, unlike the usually slightly curved paths of natural extension movements (although there are also straight paths in natural movements) (Gutfreund et al. 1998).
The tangential velocity profile was calculated using the bend point path (Fig. 7C). As in the natural movement, the profile of the simulated reaching movement can also be divided into 3 phases: 1) a minor decrease in velocity; 2) a monotonic increase in velocity; and 3) a plateau and a minor decrease in velocity at the end. Figure 8C shows velocity profiles of simulated movements using the nonlinear muscle model (the linear muscle model profiles were almost identical). Note that the velocity profiles differ in the maximal velocity and time to peak, as a function of the activation signal parameters ac_{a} and ac_{t}. The relation between ac_{a} and ac_{t} and the control of the arm extension movement are examined in detail in paper II.
COMPARISON OF THE TWO MUSCLE MODELS.
Simulating the extension movements using the 2 muscle models yielded almost identical arm shapes (Fig. 6, C and D) and velocity profiles. The only significant difference was that the linear muscle model was an order of magnitude faster to simulate than the nonlinear one. This difference in running time is mainly attributed to the adaptive timestep property of the numerical integration tool used. In the nonlinear muscle model, the numerical integrator required much smaller time steps than the linear case to achieve the same error bounds.
Sensitivity analysis
We examined to what extent the output of the model is sensitive to changes in the parameters used or to the modeling assumptions.
NUMBER OF SEGMENTS.
We simulated reaching movements using a model arm with either 20 or 40 segments. Figure 8, C and D shows the bend point velocity profiles of simulated reaching movements with constant ac_{a} of 0.6 and ac_{t} values ranging between 0.6 and 1.2 s in intervals of 0.1 s. Movements of the 20segment (Fig. 8C) and the 40segment arms (Fig. 8D) showed similar 3phased velocity profiles. Their maximal velocities were almost identical. Movements of the 40segment arm showed slightly lower minimal velocities, a small shift (of ∼0.05 s) of the velocity profiles to the right, and a small elongation of the third phase (of 0.05 to 0.1s). Simulating reaching movements of a 40segment arm took approximately twice as long as for the 20segment case. We checked other conditions (activation times and activation amplitudes) and concluded that the results presented in this study are not unique to an arm with 20 segments.
DRAG COEFFICIENTS.
These values dictate how dominant the drag forces are. Changes in the drag coefficients caused changes in the shape of the arm during movement. We simulated a reaching movement with different values of C_{per} and C_{tan} and compared the shape of the arm to the results from a standard simulation (Fig. 9, standard), i.e., the one with the experimentally derived drag coefficients. The curved surface Ψ (Fig. 9, middle) depicts the difference between the standard simulation and other simulations where the surface represents the values of where q_{i}(t) is the position of mass i in time t and the asterisk marks the positions obtained in the standard simulation.
The values of the 2 axes (Fig. 9, middle) range from 0 to 4 times the experimentally found values of C_{per} (1.013) and C_{tan} (0.0256). The difference surface Ψ(C_{per}, C_{tan}) is color coded from blue (a low difference) to red (a high difference). When both coefficients were very low, the tip of the arm overshot (Fig. 9, C_{per} = 0, C_{tan} = 0), resulting in a movement very different from the standard one (and from the natural reaching movements). A low C_{tan} had a smaller influence, causing the distal part of the arm (from bend to tip) to wiggle and the arm to be straighter toward the end of the movement (Fig. 9, C_{tan} = 0). Tip oscillations occurred only when C_{tan} was low. A low C_{per} resulted in a lower radius of curvature for the bend, a straighter arm toward the end of the movement, and a small overshoot at the end (Fig. 9, C_{per} = 0). High values of C_{tan} (Fig. 9, C_{tan} = 4) resulted in the arm being more curved toward the end of the simulation—as if the arm musculature did not have enough force to straighten the arm, as happens in normal reaching movements. High values of C_{per} (Fig. 9, C_{per} = 4) had a much smaller effect. Based on the above and the orientation of the low regions of the difference surface we conclude that the model is nearly 2.5 times more sensitive to the value of C_{tan} than to C_{per}. It is also clear that the model is not invariant to the ratio between C_{per} and C_{tan} because such invariance would be reflected in a symmetry of the difference surface along the main diagonal.
VELOCITY OF THE ACTIVATION SIGNAL.
The third phase of the velocity profiles of natural reaching movements (Fig. 8A) decreases more strongly than that of the simulated profiles (Fig. 8C). This could result from the biological activation signal differing from what we used or from different physical properties distally and proximally along the arm.
We compared the velocity profiles of simulated reaching movements using either an activation signal with constant velocity or with a bellshaped velocity profile (Fig. 5C). The simulations had constant activation amplitudes of 0.6 and activation times ranged from 0.6 to 1.3 s. The velocity profiles of the movements whose activation signal had a bellshaped velocity profile (Fig. 8E) were very different from those of the constantvelocity case (Fig. 8C). They had higher maximal velocities (resulting from higher signal velocities) and a large shift in the time to peak (resulting from the existence of a peak in the velocity of the activation signal). In many of the bellshaped velocity profiles of the activation signal, the maximal signal velocity was too high for the activation signal to successfully propel the bend and the arm was not fully straightened as it is in natural movements. Activation signals with constant velocities were able to straighten the arm in most cases, producing movements that achieved the natural range of bend velocities and fitted the observed behavior much better.
TERMINATION OF THE ACTIVATION SIGNAL.
We simulated reaching movements with an activation signal with constant velocities that gradually decreased in amplitude from maximal activation to zero, terminating at different locations along the arm in different simulations. This velocity decrease took place over approximately 0.4 of the arm length, its position varying from the midpoint of the arm to closer to the tip. The velocity profiles demonstrated a variable third phase that ranged from a large decrease in velocity (early termination) to a velocity plateau (late or no termination) (Fig. 8F).
DISTRIBUTION OF MUSCLE PROPERTIES.
We simulated reaching movements in which the ratio of the muscle forces of the longitudinal muscles to those of the transverse muscles varied from 0.5 to 2.0. Simulations with larger longitudinal muscle forces resulted in shorter (and thicker) arm shapes, whereas those with stronger transverse muscles resulted in longer (and thinner) arm shapes. The change in final arm length as a function of the ratio of muscle forces was only 20% (from 90 to 110% of the length when the ratio was 1:1) (Fig. 10A). This means that arm length in the model was not very sensitive to the ratio of the longitudinal to transverse muscles forces.
We also varied the ratio of the forces of the dorsal and ventral muscles from 0.5 to 2.0 (Fig. 10B). The final arm shape was significantly affected; stronger dorsal muscles resulted in an upward curving arm and stronger ventral muscles caused the arm to curve downward. In natural reaching movements, the arm is straightened during the movement (Fig. 6A).
DISCUSSION
The similarity between the two muscle models
The 2 muscle models used here yielded highly similar reaching movements. During the simulations, the muscle lengths usually remain between 0.7l_{0}^{m} to 1.1l_{0}^{m} (paper II, Fig. 3), where the nonlinear model force–length relationship is nearly linear (Fig. 4). Thus the (active and passive) forces produced by the nonlinear muscle model are similar to those produced by the linear muscle model. Is this what happens in the octopus arm? It is not known within what range of muscle lengths the arm operates during reaching or other movements.
The drag forces
Our experimentally obtained value for the perpendicular component of the drag force conforms to the results of an experiment on a circular cylinder described in Streeter (1962), who gives the drag coefficient of the perpendicular component versus Reynolds number. For the range of Reynolds number that characterizes the octopus' arm (≅10^{3}), the perpendicular drag coefficient C_{per} is constant with a value around 1.0, similar to our result of C_{per} = 1.013 ± 0.055. The circular cylinder is thus a good approximation for the shape of an octopus arm segment with respect to the drag forces in water.
The large difference between the normal and tangential drag coefficients may have played a major role in the evolution of a strategy for moving the flexible octopus arm from one location to another. During simulated reaching movements, the effect of the nearly 50fold larger perpendicular drag coefficient is minimized to such an extent that it approximately equals the tangential drag forces. Not all octopus arm movements reduce drag forces in this way. For example, during fetching food to the mouth, the arm forms an articulatedlike structure with parts of the arm oriented normal to the direction of movement (Sumbre et al. 2005). An octopus may also rotate an already extended arm around the arm's base to a new orientation to catch an object. However, this is very rare; octopuses usually pull back the extended arm and extend an arm again (Gutfreund et al. 1996, Fig. 1).
Comparing the simulations to natural reaching movements
The kinematic features of simulated movements were compared to the natural reaching movements, such as arm shape, bend point path, and its velocity profile. Generally, the simulation produced a bend propagation movement that straightened the arm, similar to the natural and evoked movements. However, some finer details were different. Some of these discrepancies may be used to finetune the model or to suggest different alternatives for modeling the behavior. However, we stress that the description of the natural movement is incomplete, especially regarding the neural activation command to the muscles. The nature of neural activity at the beginning and the end of a reaching movement is unknown and the measurements of activation signal velocities are incomplete. Also unknown is the extent to which the muscle groups receive different activations.
Taking the similarities and discrepancies one by one we discovered the following. 1) In natural reaching movements, the arm is almost straight toward the end of the movement. In the simulations, different activation levels of the dorsal and ventral muscles left the arm curved (Fig. 10B), whereas equal activation resulted in an almost straight arm. Our modeling thus suggests equal activation of dorsal and ventral longitudinal muscles in the natural reaching movement. 2) We used a sigmoid activation wave (a smoothed step function) that elevated muscle activity from zero to a maximal level, where it remained until the end of the simulation. This sustained muscle activity stiffens the proximal part of the arm and keeps it pointing in one direction, resisting gravity and other perturbing forces. In contrast, the EMG patterns show that the sustained muscle activity is lower than the maximal level. Using such a scheme in the model produced similar results to those produced by the sigmoid shape. The amplitude of the sustained muscle activity can thus be lower than maximal while achieving the same stability of the activated part. 3) The bend paths of simulated movements are straight and those of natural movements are usually slightly curved. This can be attributed to the rotation of the arm around its base during natural movements that usually occurs at the beginning of the movement and is variable from movement to movement (unpublished results). We did not model base rotations. 4) The first phase of the velocity profile of natural movements is variable with some profiles having a noticeable velocity before the second phase. In contrast, the velocity profiles of the simulated movements show a low velocity during the first phase. This difference is mainly explained by the use of tangential velocity profiles to describe the movement of the bend along the arm. Tangential velocity is a scalar value that measures velocity in the direction of the movement. If a bend propagates along the arm creating a straight spatial path, its tangential velocity and propagation velocity will be the same. Even if the bend location does not change at all along the arm, however, movements of the arm in space will give positive tangential velocities of the bend. Rotating the arm around its base (common in natural movements) creates a curved bend path that increases the tangential velocity beyond values deriving from bend propagation alone. In contrast, the bend paths of simulated movements are straight and their tangential velocities are directly related to the bend propagations. 5) The third phase of natural velocity profiles, usually showing a velocity decrease, is much more variable than in the simulations, where movements with a constant velocity activation signal show a plateau followed by a small decrease. The difference could result from variability in the neural command, such as the timing of signal termination. Termination of the activation signal causes a decrease in muscle activity and any decrease in the muscle forces propelling the bend can cause the observed decrease in bend velocity (Fig. 8F). Anatomically, the muscle cells have similar properties along the arm (unpublished observation). These data do not support the decrease in velocity in natural movements as attributed to weaker muscles around at the tip of the arm.
Based on the above, we hypothesize that the velocity of natural activation signals probably lies somewhere between the constant and bellshaped velocity forms we examined and that the amplitude of the activation command decreases toward the end of the movement. The large variability in natural movements does not allow us to draw a more precise conclusion.
Finally, this paper describes a dynamic model of the octopus arm that is general enough for modeling other movements in the octopus, as well as other in muscular hydrostats. Exploring the hypothesis of a “stiffening wave,” we have shown that an extension movement can be achieved by a very simple control scheme that minimizes the computational complexity of motor control (which is an important factor for a hyperredundant arm). The stiffening wave was also found to be an efficient way of moving in water. From here, we can speculate on the evolutionary pressures that led to the development of the octopus' neuromuscular system.
In paper II, we use the model developed here to study the biomechanical details of the extension movements and its control. We show why equal activation of the muscles of an arbitrarily shaped segment transforms it into a symmetric rectangular shape. Such shape changes in several consecutive segments result in the arm straightening. We then focus on the properties of the activation signal that creates a wave of muscle contraction. We show that the 2 properties of the activation signal, its amplitude and its traveling time, are related and the motor control system may use this relation to choose the minimal forces needed to extend an arm according to a desired kinematic plan. We also show that it is possible to use higher forces to resist external perturbations while keeping the same kinematics. Therefore, our model suggests that the octopus can independently control 2 major aspects of arm extension movements: the kinematics of the movement and its resistance to perturbations. We do not know whether the behaving octopus actually performs according to this prediction.
APPENDIX A
Matrices and vectors definitions
The model arm is composed of n mass pairs (total of 2n masses), creating n − 1 segments. M is a diagonal mass matrix where each mass is represented twice (dimensions: 4n × 4n) where the superscript v denotes ventral and d denotes dorsal (see Fig. 1).
q is the position vector, composed of the coordinates of the 2n masses (dimensions: 4n × 1) again v denotes ventral and d denotes dorsal.
The same inner organization appears in the force vectors:
f^{m} denotes the muscle forces (4n × 1)
f^{g} denotes the combined gravity and buoyancy forces (4n × 1)
f^{w} denotes the water drag forces (4n × 1)
f^{c} denotes the forces arising from the fixed volume constraints (4n × 1)
s is the vector of areas of the n − 1 segments
The area of each segment was calculated by decomposing the quadrilateral shape into 2 general triangles and finding their areas
s^{c} is the vector of values of the constant areas for the n − 1 segments
D is a matrix (n − 1 × 4n) transforming the position vector q to the vector s, which are the areas of the n − 1 segments of the arm
APPENDIX B
Finding the matrix C in the constraint force by the Lagrange formulation
The constraint force should counteract all changes in shape that lead to changes in the area of a segment. We determine the constraint force by finding the pressures in each compartment and using the Lagrange formulation with Lagrange multipliers corresponding to the pressures:
The Lagrangian is^{9} where the E_{k} is the kinetic energy where s_{i} is the area of segment i (see the definition of the vector s above), s_{i}^{c} is the constraint on the area, λ_{i} is the Lagrange multiplier (corresponding to the pressure of compartment i), q is the vector (x, y) of generalized coordinate, and q̇ is the time derivative of q.
The Euler–Lagrange equation is (B1) By obtaining the derivatives, we obtain explicit expressions for the 2 terms of Eq. B1. For example, for the ith ventral mass we have Note that s_{i}_{−1} and s_{i} are the only segments whose areas depend on x_{i}^{v}.
Substituting into Eq. B1 we get where m_{i}^{v}ẍ_{i}^{v} is exactly the constraint force for the x_{i}^{v} coordinate. The solutions for the other coordinates (x_{i}^{d}, y_{i}^{v}, and y_{i}^{d}) are obtained in the same way.
Writing the constraint force vectorially, we get f^{c} (4n × 1) Let us write f^{c} as Cλ, where λ = (λ_{1} · · · λ_{n}_{−1}). Note that λ equals the pressure vector p from Eq. 3.
The matrix C (4n × n − 1) is
Achieving constraint equations at the acceleration level
Here we use index notation^{10} for its brevity.
The constraint equation is where s^{c} is the vector of constant areas of the segments.
Differentiating the constraint equation with respect to time, we have Differentiating again with respect to time Using the chain rule we have of which the first term equals
Every element in D depends linearly on q, so ∂D_{ik}/∂q_{l} is a constant element and differentiating it with respect to time leads to 0, leaving only the last term. By rearranging the equation we have Defining (45) and we arrive at (B2)
Finding the compartmental pressures
Multiplying Eq. 1 by M^{−1} we get (B3) Writing f^{c} as Cp (Eq. 3) and multiplying Eq. B3 by G results in (B4) We then substitute Gq̈ from Eq. B2, multiply Eq. B4 by (GM^{−1}C)^{−1} and rearrange (B5) leading to f^{c}.
The equations of motion can now be solved using numerical integration (we use Matlab's ode23, a mediumorder variable step, Runge–Kutta method). The 4 forces and the masses are known, so the acceleration (q̈) can be found directly from Newton's second law. Integrating the accelerations leads to the velocities (q̇), and a second integration to the positions (q), which are used in the next time step to calculate the new forces.
GRANTS
This work was supported by Defense Advanced Research Projects Agency Grant N6600103R8043, Israel Science Foundation Grant 580/02, and the Moross laboratory. T. Flash is an incumbent of the Dr. Hymie Moross Professorial chair.
Acknowledgments
We thank Dr. Jenny Kien for suggestions and editorial assistance; I. Givoni, A. Shabtay, and R. Adam for valuable discussions; and G. Sumbre for advice and assistance.
Footnotes
↵* Y. Yekutieli and R. SagivZohar contributed equally to this work.
↵1 Reynolds number (Re) is the dimensionless ratio of inertial to viscous forces of bodies moving in a fluid. Re is given by a characteristic velocity (v) of the moving body, multiplied by the characteristic length (l) of the moving body and divided by υ, the kinematic viscosity of the fluid. The Reynolds number indicates flow regime (laminar or turbulent flow) [see Streeter (1962), p. 161–168, and Vogel (1981), p. 70, 135].
↵2 Dimensionless analysis is a technique used in the physical sciences and engineering to reduce physical properties such as acceleration, viscosity, energy, and others to their fundamental dimensions of length (L), mass (M), and time (T). For example, a force dimension is [F] = [ML/T^{2}] (see Batchelor 1967; Streeter 1962).
↵3 Steady flow occurs when conditions at any point in the fluid do not change with time. The definition of a steady flow is generalized (using averaged velocity, pressure, and temperature) to account for the fluctuations in turbulent flow (see Streeter 1962).
↵4 In a crossflow experiment, the body is static and the direction of flow is perpendicular to the long axis of the body.
↵5 In a streamwiseflow experiment, the body is static and the direction of flow is tangential to the long axis of the body.
↵6 Relaxed length is the length of a muscle that is not activated and is not subjected to external forces. A muscle set to this length can produce force only if it is activated. For muscle lengths longer than the relaxed length, a resistive force can be detected without activation (Fig. 6A).
↵7 Minimum jerk polynomials were used by Flash and Hogan (1985) to quantify the smoothness of pointtopoint reaching movements in humans. Here we use the following polynomial to describe the velocity profile of the activation signal where L is the total distance that the signal travels (arm length) and T is the total traveling time.
↵8 Tangential velocity is the velocity in the direction of movement. It is calculated from the trajectory
↵9 The potential energy of gravity and buoyancy was ignored because it amounts to a small, timeindependent force directed along the yaxis, which is practically identical for the masses of a segment.
↵10 A vector has one subscript (q_{k}) and a matrix has two (D_{ik}). Summation is denoted by a repeated index: D_{ik}q_{k} = s_{i}^{c} means ∑_{k} D_{ik}q_{k} = s_{i}^{c}, which is the ith element of the matrix equation Dq = s_{i}^{c}.
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
 Copyright © 2005 by the American Physiological Society