| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
1Department of Neurobiology and 2Interdisciplinary Center for Neural Computation, Hebrew University, Jerusalem; and 3Department of Computer Science and Applied Mathematics, Weizmann Institute of Science, Rehovot, Israel
Submitted 6 July 2004; accepted in final form 2 March 2005
 |
ABSTRACT |
|---|
|
|
|
INTRODUCTION |
|---|
|
). 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.
|
; 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.
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.) 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. 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. ; Wadepuhl and Beyn 1989). 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 MECHANICAL MODEL—AN OVERVIEW |
|---|
|
for the arrangement of the musculature of the octopus arm.]
|
|
2 or 3 of a physical quantity like velocity
m, m > 3
, 
, respectively
Model assumptions and limitations
|
|
THE MODEL IN DETAIL |
|---|
|
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:
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) |
). 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) |
is density, 1,022 kg/m3 for seawater and 1,042 kg/m3 for an octopus arm; Varm is a 4n x 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
 | (3) |
(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
 | (4) |
.)
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
 | (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)
 | (6) |
 | (7) |
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
 |
 |
).
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
 |
 |
, 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 x 35 x 40-cm (length x width x 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.
|
where 
are force vectors). In our experiment, the relevant forces acting on the segment were confined to the vertical direction.
 | (8) |
 | (9) |
is Archimedes' force acting on the segment and fAlo = Vlo 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 2
l, where is the mean segment radius of the arm and l is the segment length. The sectional surface area Sa is 2
l.
For the perpendicular component Re = 2
per/. For the tangential component Re = ltan/, where l is the segment length, is the segment average velocity, is the segment average radius, and is the kinematic viscosity ( = 1 x 10–6 m2/s).
ANALYSIS OF THE EXPERIMENT. 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.
|
Pacper (Fig. 3A and Eq. 6), and for the tangential component the slope of the linear fit is 0.5Sactan (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 cper (from Eq. 6)
 |
For the tangential drag coefficient ctan (from Eq. 7)
 |
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.
and Zajac (1989). It includes nonlinear force–length and force–velocity relations. ). 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) |
[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.
|
CHOICE OF PARAMETERS IN THE NONLINEAR MUSCLE MODEL.
The relaxed cross-sectional area differs between the longitudinal and transverse muscles (Fig. 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).
|
, 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).
|
; Milligan et al. 1997). Peak stress in a brief isometric tetanus was 1.31 x 105 N/m2 in squid tentacle muscles and 4.68 x 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 x 105 N/m2 (Lowy 1962). Our measurements on a lateral group of longitudinal octopus arm muscles gave 4.0 x 104 N/m2 (H. Matzner and B. Hochner, unpublished results). The maximal active muscle stress used in the model here is 8.0 x 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
 |
). 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
 | (11) |
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. PARAMETER CHOICE FOR THE LINEAR MUSCLE MODEL. 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 x 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 x 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
 |
 |
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. 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 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
 | (12) |
(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)
 | (13) |
). 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).
|
|
RESULTS |
|---|
|
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
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 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).
|
. 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.
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 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.
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 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.
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 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
 |
|
(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. 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 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.
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 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.
|
|
APPENDIX A |
|---|
|
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 x 4n)
 |
q is the position vector, composed of the coordinates of the 2n masses (dimensions: 4n x 1)
 |
The same inner organization appears in the force vectors:
The area of each segment was calculated by decomposing the quadrilateral shape into 2 general triangles and finding their areas
 |
|
|
APPENDIX B |
|---|
|
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
 |
 |
 |
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
 | (B1) |
 |
 |
Substituting into Eq. B1 we get
 |
iv 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 x 1)
 |
, where = (1 · · · n–1). Note that equals the pressure vector p from Eq. 3.
The matrix C (4n x n – 1) is
 |
Achieving constraint equations at the acceleration level
Here we use index notation10 for its brevity.
The constraint equation is
 |
Differentiating the constraint equation with respect to time, we have
 |
 |
 |
 |
 |
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
 |
 | (45) |
 |
 | (B2) |
Finding the compartmental pressures
Multiplying Eq. 1 by M–1 we get
 | (B3) |
 | (B4) |
from Eq. B2, multiply Eq. B4 by (GM–1C)–1 and rearrange
 | (B5) |
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.
|
|
GRANTS |
|---|
|
|
|
ACKNOWLEDGMENTS |
|---|
|
|
|
FOOTNOTES |
|---|
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.
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
 |
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, 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.
Address for reprint requests and other correspondence: T. Flash, The Faculty of Mathematics and Computer Science, The Weizmann Institute of Science, POB 26, Rehovot 76100, Israel (E-mail tamar.flash{at}weizmann.ac.il)
|
|
REFERENCES |
|---|
|
Alscher C and Beyn WJ. Simulating the motion of the leech: a biomechanical application of DAEs. Numer Algorithms 19: 1–12, 1998.[CrossRef]
Bogacki P and Shampine LF. A 3(2) pair of Runge–Kutta formulas. Appl Math Lett 2: 1–9, 1989.
Chiel HJ, Crago PE, Mansour JM, and Hathi K. Biomechanics of a muscular hydrostat: a model of lapping by a reptilian tongue. Biol Cybern 67: 403–415, 1992.[CrossRef]
Chirikjian GS and Burdick JW. A modal approach to hyper-redundant manipulator kinematics. IEEE Trans Robot Autom 10: 343–354, 1994.[CrossRef]
Curtin NA, Woledge RC, and Bone Q. Energy storage by passive elastic structures in the mantle of sepia officinalis. J Exp Biol 203: 869–878, 2000.[Abstract]
Drushel RF, Neustadter DM, Hurwitz I, Crago PE, and Chiel HJ. Kinematic models of the buccal mass of Aplysia californica. J Exp Biol 201: 1563–1583, 1998.
Drushel RF, Sutton GP, Neustadter DM, Mangan EV, Adams BW, Crago PE, and Chiel HJ. Radula-centric and odontophore-centric kinematic models of swallowing in Aplysia californica. J Exp Biol 205: 2029–2051, 2002.
Ekeberg Ö. A combined neuronal and mechanical model of fish swimming. Biol Cybern 69: 363–374, 1993.[CrossRef]
Flash T and Hogan N. The coordination of arm movements: an experimentally confirmed mathematical model. J Neurosci 5: 1688–1703, 1985.[Abstract]
Gutfreund Y, Flash T, Fiorito G, and Hochner B. Patterns of arm muscle activation involved in octopus reaching movements. J Neurosci 18: 5976–5987, 1998.
Gutfreund Y, Flash T, Yarom Y, Fiorito G, Segev I, and Hochner B. Organization of octopus arm movements: a model system for studying the control of flexible arms. J Neurosci 16: 7297–7307, 1996.
Jordan CE. A model of rapid-start swimming at intermediate Reynolds number: undulatory locomotion in the chaetognath Sagitta elegans. J Exp Biol 163: 119–137, 1992.
Jordan CE. Coupling internal and external mechanics to predict swimming behavior: a general approach. Am Zool 36: 710–722, 1996.
Kier WM. The arrangement and function of molluscan muscle. In: The Mollusca, Form and Function, edited by Trueman ER, Clarke MR, and Wilbur KM (editor-in-chief). New York: Academic Press, 1988, vol. 11, p.211–252.
Kier WM. Hydrostatic skeletons and muscular hydrostats. In: Biomechanics (Structures and Systems): A Practical Approach, edited by Biewener AA. Oxford, UK: IRL Press, 1992, p. 205–231.
Kier WM and Curtin NA. Fast muscle in squid (Loligo pealei): contractile properties of a specialized muscle fiber type. J Exp Biol 205: 1907–1916, 2002.
Kier WM and Smith KK. Tongues, tentacles and trunks: the biomechanics of movement in muscular-hydrostats. Zool J Linn Soc 83: 307–324, 1985.
Kristan WB, Skalak R, Wilson RJA, Skierczynski BA, Murray JA, Eisenhart FJ, and Cacciatore TW. Biomechanics of hydroskeletons: lessons learned from studies of crawling in the medicinal leech. In: Biomechanics and Neural Control of Posture and Movement, edited by Winters J and Crago P. New York: Springer-Verlag, 2000, p. 206–218.
Lowy J and Millman BM. Mechanical properties of smooth muscles of cephalopod molluscs. J Physiol 160: 353–363, 1962.
Matzner H, Gutfreund Y, and Hochner B. Neuromuscular system of the flexible arm of the octopus: physiological characterization. J Neurophysiol 83: 1315–1328, 2000.
Milligan BJ, Curtin NA, and Bone Q. Contractile properties of obliquely striated muscle from the mantle of squid (Alloteuthis subulata) and cuttlefish (Sepia officinalis). J Exp Biol 200: 2425–2436, 1997.[Abstract]
Neustadter DM, Drushel RF, Crago PE, Benjamin WA, and Chiel HJ. A kinematic model of swallowing in Aplysia californica based on radula/odontophore kinematics and in vivo magnetic resonance images. J Exp Biol 205: 3177–3206, 2002.
Skierczynski BA, Wilson RJA, Kristan WB, and Skalak R. A model of the hydrostatic skeleton of the leech. J Theor Biol 181: 329–342, 1996.[CrossRef][Web of Science][Medline]
Streeter VL. Fluid Mechanics (3rd ed.). Kogakusha, Japan: McGraw-Hill, 1962, chapt. 3.
Sumbre G, Fiorito G, Flash T, and Hochner B. Precise motor control of flexible arms. Nature 433: 595–596, 2005.[CrossRef][Medline]
Sumbre G, Gutfreund Y, Fiorito G, Flash T, and Hochner B. Control of octopus arm extension by a peripheral motor program. Science 293: 1845–1848, 2001.
Takanashi N, Choset H, and Burdick JW. Experimental results of sensor based planning for hyper-redundant manipulators. Proc IEEE/IROS Yokohoma, Japan, 1993.
Van Leeuwen JL and Kier WM. Functional design of tentacles in squid: linking sarcomere ultrastructure to gross morphological dynamics. Philos Trans R Soc Lond B Biol Sci 352: 551–571, 1997.[CrossRef]
Vogel S. Life in Moving Fluid: The Physical Biology of Flow Boston, MA: Willard Grant Press, 1981, p. 70, 135.
Wadepuhl M and Beyn WJ. Computer simulation of the hydrostatic skeleton. the physical equivalent, mathematics and application to worm-like forms. J Theor Biol 136: 379–402, 1989.[CrossRef][Web of Science][Medline]
Walker ID. Some issues in creating "invertebrate" robots. Proc Int Symp Adaptive Motion of Animals and Machines, Montreal, Canada, 2000.
White FM. Viscous Fluid Flow. New York: McGraw-Hill, 1974.
Wilson JF, Mahajan U, Wainwright SA, and Croner LJ. A continuum model of elephant trunks. J Biomech Eng 113: 79–84, 1991.[Web of Science][Medline]
Yekutieli Y, Sagiv-Zohar R, Hochner B, and Flash T. Dynamic model of the octopus arm. II. Control of reaching movements. J Neurophysiol 94: 1459–1468, 2005.
Yu SN, Crago PE, and Chiel HJ. Biomechanical properties and a kinetic simulation model of the smooth muscle I2 in the buccal mass of Aplysia. Biol Cybern 81: 505–513, 1999.
Zajac FE. Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. CRC Crit Rev Biomed Eng 17: 359–411, 1989.[Web of Science][Medline]
This article has been cited by other articles:
|
|
 |
|
  R. J. Gilbert, V. J. Napadow, T. A. Gaige, and V. J. Wedeen Anatomical basis of lingual hydrostatic deformation J. Exp. Biol., December 1, 2007; 210(23): 4069 - 4082. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 Y. Yekutieli, R. Mitelman, B. Hochner, and T. Flash Analyzing Octopus Movements Using Three-Dimensional Reconstruction J Neurophysiol, September 1, 2007; 98(3): 1775 - 1790. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
C. L. Huffard Locomotion by Abdopus aculeatus (Cephalopoda: Octopodidae): walking the line between primary and secondary defenses J. Exp. Biol., October 1, 2006; 209(19): 3697 - 3707. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
L. Blackburn AGILE ANIMALS J. Exp. Biol., November 1, 2005; 208(21): iv - iv. [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Visit Other APS Journals Online |
TOP
ABSTRACT
INTRODUCTION
