Dynamic Model of the Octopus Arm. I. Biomechanics of the Octopus Reaching Movement

Yoram Yekutieli, Roni Sagiv-Zohar, Ranit Aharonov, Yaakov Engel, Binyamin Hochner, Tamar Flash


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.


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 fluid-filled 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.

  1. The models are static or quasi-static, i.e., they do not describe the full dynamics of motion: Wadepuhl and Beyn (1989) developed a quasi-static model of wormlike forms; Skierczynski et al. (1996) and Kristan et al. (2000) developed a quasi-static 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.)

  2. 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.

  3. 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.

  4. 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).

  5. 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 2-dimensional (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 z-part 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 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.]

FIG. 1.

2n masses are arranged in n pairs, each consisting of one ventral and one dorsal mass and are connected by longitudinal and transverse muscles. Alongitudinal and Atransverse are the relaxed cross-sectional areas that are perpendicular to the longitudinal and transverse muscles, respectively.

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:

a scalar
a vector in ℜ2 or ℜ3 of a physical quantity like velocity
a vector in ℜm, m > 3
a scalar average value
a matrix
xi, xij
elements of a vector and matrix, respectively
time derivatives of x and , respectively

Model assumptions and limitations

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

  2. 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).

  3. 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.


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:

  1. The internal forces generated by the arm muscles (fm).

  2. The vertical forces resulting from the combined influences of gravity and arm buoyancy in seawater (Archimedes' law of buoyancy) (fg).

  3. The drag produced by arm motion through the water (fw).

  4. The internal forces responsible for maintaining the constant volume constraints (fc).

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 Math(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 Math(2) where ρ is density, 1,022 kg/m3 for seawater and 1,042 kg/m3 for an octopus arm; Varm is a 4n × 4n matrix of segment volumes; and g is the gravitational acceleration of the Earth.

The total weight fg 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 fc is linear and can be written compactly as Math(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 fc cannot be directly determined. The constraint force fc must be derived indirectly from the constant volume constraints (see appendix B). The result of the derivation is the pressure vector p Math(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 fc (Eq. 3). All forces are now expressed explicitly and the equations of motion can be solved Math(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: vper in the direction perpendicular to the long axis of the segment and vtan in the direction tangential to the long axis. Similarly, the water drag forces acting on the segment can be decomposed into a perpendicular component dper and a tangential component dtan.

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 flow3are (Vogel 1981) Math(6) Math(7) where Pa is the projected area of the segment in the plane perpendicular to the direction of ν⃗per and Sa 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 cper and ctan express the perpendicular and tangential drag coefficients, respectively, which are functions of Reynolds number as follows Math Math The functions f and g are of the form u1 + u2(Re)u3, where u1, u2, and u3 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 (cper) from curves fitted for cylinders in cross flow, 4 and a tangential drag coefficient (ctan) 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 Math Math where Reper, 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 Retan 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, dper and dtan. First we calculated the real projected area Pa and sectional area Sa. We then used Eq. 6 and 7 to obtain the experimental drag coefficients cper and ctan. 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 × 40-cm (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 S-VHS video camera with 40 ms between adjacent frames.

FIG. 2.

An amputated octopus arm was dropped in seawater under different gravitational loads and its falls were used to calculate the 2 components of the drag force; the segment accelerated and reached a constant velocity, from which the drag forces could be calculated using the force balance equation. A: snapshots of horizontal falls; the experiment aimed to measure the velocity component perpendicular to the long axis of the arm segment. B: snapshots of vertical falls to find the tangential velocity, i.e., parallel to the long axis of the arm segment. Time is shown at the top of each frame. White arrows mark direction of the different forces acting on the arm segment: fAse,fAlo, buoyancy and gravitation of the segment and of the load, respectively; dper, dtan, perpendicular and tangential drag forces, respectively; Wse, Wlo, weight of the segment and of the load, respectively.


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 Math 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 Math(8) Math(9) where fAse = Vseρg⃗ is Archimedes' force acting on the segment and fAlo = Vloρg⃗ is Archimedes' force acting on the load. The volumes of the segment and the load are Vse and Vlo, respectively, and the weight of the segment and the load are Wse and Wlo.

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 Pa is 2r̄l, where is the mean segment radius of the arm and l is the segment length. The sectional surface area Sa 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, is the segment average velocity, is the segment average radius, and υ is the kinematic viscosity (υ = 1 × 10−6 m2/s).


For each experiment, we plotted velocity versus time, found the leveling-off 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.

FIG. 3.

A linear fit was used to calculate the drag coefficients: A: measured drag force versus squared velocity in the vertical fall; slope is 0.41 ± 0.02. B: measured drag force vs. squared velocity in the horizontal fall; slope is 0.03 ± 0.002.

For the perpendicular component, the slope of the linear fit is 0.5ρPacper (Fig. 3A and Eq. 6), and for the tangential component the slope of the linear fit is 0.5ρSactan (Fig. 3B and Eq. 7).


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 cper (from Eq. 6) Math

For the tangential drag coefficient ctan (from Eq. 7) Math

Note that the perpendicular drag coefficient is nearly 50-fold 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.

  1. 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.

  2. 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 Math(10) where A is the cross-sectional area of the muscle, a ∈ [0, 1] is a dimensionless activation function, l(t) is muscle length, and l0m is the specific length of the muscle at which active muscle force achieves its peak value (terminology after Zajac 1989). The active stress Fa equals the normalized active force (Fig. 4A) multiplied by the maximum stress. The dimensionless normalized force–velocity relation is expressed by Fv, which is a function of v/vmax, where v = dl(t)/dt, the rate of change of muscle length with time, and vmax is the absolute value of the maximum shortening velocity. The passive stress Fp is equal to the normalized passive force (Fig. 4A) multiplied by the maximum stress.

FIG. 4.

A: passive and active force–length relationships in the nonlinear muscle model. Force is expressed relative to the peak active force and length relative to the optimal length of active force production l0m. Inactive muscle starts to resist pulling for muscle lengths larger than the relaxed length. B: nonlinear muscle model force–velocity relationship. Again, force is normalized to peak active force and velocity to the absolute value of the maximum shortening velocity vmax. C: passive and active force–length relationships in the linear muscle model, normalized as in A. Linear relations (straight lines) are superimposed on the nonlinear curves. See text for explanation.

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 Fp (Fig. 4A, “passive” curve). In a fully activated isometric contraction the muscle force is the sum of the active and passive contributions Fa and Fp (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 Fv (see Fig. 4B and explanation below).


The relaxed cross-sectional area differs between the longitudinal and transverse muscles (Fig. 1) Math Math where wr and lr are the relaxed lengths6 of the transverse and longitudinal segmental muscles, respectively; and d is the depth of the segment, which equals wr at the initial relaxed shape of the arm. The factor of in Alongitudinal 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 Fa, Fp, and l0m: 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).

FIG. 5.

Activation function used in the model. A: sigmoid activation function vs. time approximated the EMG measured by Gutfreund et al. (1998). B: activation along the arm at a specific time step for a signal propagating at a constant velocity. All muscles to the left of the top circle are fully activated. All muscles to the right of the bottom circle are not activated at all, and the muscles in between are partially activated according to the sigmoid shape. C: 2 velocity types used for the activation signal, a constant velocity (horizontal line) and a bell-shaped velocity profile. Activation signal takes the same time to move from the base to the tip of the arm in both cases, and thus both have the same averaged velocity.

The peak stress in a brief isometric tetanus was found to be 2.62 × 105 N/m2 in squid mantle muscle and 2.26 × 105 N/m2 in cuttlefish (Curtin et al. 2000; Milligan et al. 1997). Peak stress in a brief isometric tetanus was 1.31 × 105 N/m2 in squid tentacle muscles and 4.68 × 105 N/m2 in the squid arm (Kier and Curtin 2002). The maximum value observed from twitches of the funnel retractor muscles of octopus was 5.0 × 105 N/m2 (Lowy 1962). Our measurements on a lateral group of longitudinal octopus arm muscles gave 4.0 × 104 N/m2 (H. Matzner and B. Hochner, unpublished results). The maximal active muscle stress used in the model here is 8.0 × 104 N/m2.

Both the active and passive force–length curves are normalized to l0m, the length at which active muscle force is maximal. The relaxed length of the group of muscles we measured was 0.8l0m (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 Fv at the boundaries are the same as those in the relation described by Zajac (1989). We use the following relation Math where v* = v/vmax is the normalized velocity and vmax 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.5l0m s−1 for squid tentacle fibers and 1.5l0m s−1 for squid arm muscle fibers.

The linear muscle model

Every linear muscle in the model exerts the following force Math(11) where lrest 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 k0, and the maximal active spring constant of the muscle by kmax, both have dimensions of N/m.


Values of kmax, k0, and lrest: 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 lrest is the length at which both models produce no force. Its value is 0.4l0m (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 Ttotal is the sum of the passive stiffness Tpassive and the active stiffness Tactive. It was taken as the slope of the line connecting the point (lrest, 0) and the maximal force point (1, 1) in Fig. 4C. This gave a value of Ttotal as 1.34 × 105 N/m2.

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 Tpassive was found to be 2,000 N/m2. The active stiffness Tactive is the difference between Ttotal and Tpassive and its value was 1.32 × 105 N/m2. The maximal spring constant kmax equals Tactive times the muscle's relaxed cross-sectional area, divided by its relaxed length. Thus, in Eq. 11 the calculation of kmax is different for a longitudinal and a transverse muscle.

For longitudinal muscles Math whereas for the transverse muscles Math The passive component of the muscle force (k0 in Eq. 11) is similarly computed using the value of Tpassive instead of Tactive. The property α of Eq. 11 is calculated analogously to kmax, where Tactive is substituted by α0 = 9 Ns/m2, 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.


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 bell-shaped 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 Math(12) where aca 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 i0 is the phase-shift parameter that is equal for all segments and was set at i0 = 1. Once the activation function f is chosen, the only remaining control parameters are aca, i0, and τ.

One such activation function is an inverted sigmoid (Fig. 5) Math(13) This sigmoid shape of the activation-time 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 act. 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 one-to-one relation between excitation and contraction, this delay may be absorbed into the phase-shift parameter i0.

We used a minimum jerk polynomial7 for the bell-shaped 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).


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 aca and the activation traveling time act, while keeping all other parameters fixed. Unless otherwise stated, we use the nonlinear muscle model.

Producing a reaching motion


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).

FIG. 6.

A: sequence of frames showing a freely behaving octopus reaching toward a target. B: electrically evoked bend propagation in the denervated arm of a decerebrated animal. C: sequence from a simulated arm extension using the linear muscle model. Red represents fully activated arm segments. Blue represents all segments not yet activated. D: as in C using the nonlinear muscle model.


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.

FIG. 7.

Simulated arm extension with an activation signal taking 1.0 s to travel along the whole arm. A: 6 traces of arm position from the simulation, showing changes over time. Red represents fully activated arm segments; blue represents all segments not yet activated. Bend point position (arrows) was marked manually. To estimate the variance in the marking process, we marked the bend point position of the same simulation 9 times. B: path of the bend point is almost straight (plotted as mean ± 1SD). C: Tangential velocity profile of the bend point (plotted as mean tangential velocity profile ± 1SD).

Bend point path and velocity profile


Octopus reaching movements are characterized by an invariant tangential velocity8 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).

FIG. 8.

A and B: velocity profiles of natural and evoked movements in denervated arms C: velocity profiles of the simulated reaching movements of a 20-segment arm using a constant velocity for the activation signal. Activation traveling time (act) ranges from 0.6 s (top curve) to 1.2 s (bottom curve). All profiles have the same activation amplitude (aca) of 0.625. D: as in C but using a 40-segment arm. E: as in C but using an activation signal with a bell-shaped velocity profile. F: as in C but all movement had the same act (0.58 s) and the activation signal terminated at different locations along the arm, from half the length of the arm (bottom velocity profiles) to the arm tip (top velocity profiles).


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 aca and act. The relation between aca and act and the control of the arm extension movement are examined in detail in paper II.


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 time-step 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.


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 aca of 0.6 and act values ranging between 0.6 and 1.2 s in intervals of 0.1 s. Movements of the 20-segment (Fig. 8C) and the 40-segment arms (Fig. 8D) showed similar 3-phased velocity profiles. Their maximal velocities were almost identical. Movements of the 40-segment 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 40-segment arm took approximately twice as long as for the 20-segment 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.


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 Cper and Ctan 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 Math where qi(t) is the position of mass i in time t and the asterisk marks the positions obtained in the standard simulation.

FIG. 9.

Effect of the values of the 2 drag coefficients on the reaching movement is demonstrated using a difference measure (colored surface) as a function of Cper and Ctan. Each point on the surface represents the difference between a simulation that uses these specific drag coefficients and the standard simulation (the one with the experimentally found drag coefficients). Evolution of arm shape during movement is shown for 7 combinations of drag coefficients.

The values of the 2 axes (Fig. 9, middle) range from 0 to 4 times the experimentally found values of Cper (1.013) and Ctan (0.0256). The difference surface Ψ(Cper, Ctan) 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, Cper = 0, Ctan = 0), resulting in a movement very different from the standard one (and from the natural reaching movements). A low Ctan 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, Ctan = 0). Tip oscillations occurred only when Ctan was low. A low Cper 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, Cper = 0). High values of Ctan (Fig. 9, Ctan = 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 Cper (Fig. 9, Cper = 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 Ctan than to Cper. It is also clear that the model is not invariant to the ratio between Cper and Ctan because such invariance would be reflected in a symmetry of the difference surface along the main diagonal.


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 bell-shaped 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 bell-shaped velocity profile (Fig. 8E) were very different from those of the constant-velocity 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 bell-shaped 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.


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).


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.

FIG. 10.

A: changes of final arm length as a function of the ratio of muscle force between the longitudinal and the transverse muscles. Length is expressed as percentage of the length at 1:1 force ratio between the 2 muscle groups. B: final arm shape as a function of the indicated ratio of dorsal to ventral muscle force.

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).


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.7l0m to 1.1l0m (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 (≅103), the perpendicular drag coefficient Cper is constant with a value around 1.0, similar to our result of Cper = 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 50-fold 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 articulated-like 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 fine-tune 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 bell-shaped 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.


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) Math 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) Math again v denotes ventral and d denotes dorsal.

The same inner organization appears in the force vectors:

  • fm denotes the muscle forces (4n × 1)

  • fg denotes the combined gravity and buoyancy forces (4n × 1)

  • fw denotes the water drag forces (4n × 1)

  • fc 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 Math

  • sc 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


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 Math is9 Math where the Ek is the kinetic energy Math where si is the area of segment i (see the definition of the vector s above), sic 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 is the time derivative of q.

The Euler–Lagrange equation is Math(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 Math Math Note that si−1 and si are the only segments whose areas depend on xiv.

Substituting into Eq. B1 we get Math where miviv is exactly the constraint force for the xiv coordinate. The solutions for the other coordinates (xid, yiv, and yid) are obtained in the same way.

Writing the constraint force vectorially, we get fc (4n × 1) Math Let us write fc as Cλ, where λ = (λ1 · · · λn−1). Note that λ equals the pressure vector p from Eq. 3.

The matrix C (4n × n − 1) is Math

Achieving constraint equations at the acceleration level

Here we use index notation10 for its brevity.

The constraint equation is Math where sc is the vector of constant areas of the segments.

Differentiating the constraint equation with respect to time, we have Math Differentiating again with respect to time Math Using the chain rule Math we have Math of which the first term equals Math

Every element in D depends linearly on q, so ∂Dik/∂ql 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 Math Defining Math(45) and Math we arrive at Math(B2)

Finding the compartmental pressures

Multiplying Eq. 1 by M−1 we get Math(B3) Writing fc as Cp (Eq. 3) and multiplying Eq. B3 by G results in Math(B4) We then substitute Gq̈ from Eq. B2, multiply Eq. B4 by (GM−1C)−1 and rearrange Math(B5) leading to fc.

The equations of motion can now be solved using numerical integration (we use Matlab's ode23, a medium-order variable step, Runge–Kutta method). The 4 forces and the masses are known, so the acceleration () can be found directly from Newton's second law. Integrating the accelerations leads to the velocities (), and a second integration to the positions (q), which are used in the next time step to calculate the new forces.


This work was supported by Defense Advanced Research Projects Agency Grant N66001-03-R-8043, Israel Science Foundation Grant 580/02, and the Moross laboratory. T. Flash is an incumbent of the Dr. Hymie Moross Professorial chair.


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.


  • * Y. Yekutieli and R. Sagiv-Zohar 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/T2] (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 cross-flow experiment, the body is static and the direction of flow is perpendicular to the long axis of the body.

  • 5 In a streamwise-flow 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 point-to-point reaching movements in humans. Here we use the following polynomial to describe the velocity profile of the activation signal Math 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 Math

  • 9 The potential energy of gravity and buoyancy was ignored because it amounts to a small, time-independent force directed along the y-axis, which is practically identical for the masses of a segment.

  • 10 A vector has one subscript (qk) and a matrix has two (Dik). Summation is denoted by a repeated index: Dikqk = sic means k Dikqk = sic, which is the ith element of the matrix equation Dq = sic.

  • 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.


View Abstract