Transcript A model for anguilliform swimming, or, a rod with a mind of its own

Nonlinear muscles,
viscoelasticity and body taper
in the creation of curvature waves
SIAM PDEs December 10, 2007
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Tyler McMillen
California State University, Fullerton
In collaboration with Thelma Williams and Philip Holmes.
Relative timing of activation and
movement
Curvature travels slower than activation.
Figure from Williams, et. al., J. Exp. Biol. (1989)
How are waves of curvature created and
propagated?
Why does curvature travel slower than
activation?
Outline
•
•
•
•
Elastic rod model
Resistive fluid forces
Discretization -- chain of rigid links
Muscle forces
Dynamics of an actuated elastic rod in the plane
g Wy
Geometry and inextensibility:
Momentum balances:
f Wx
f , g - contact forces (maintain inextensibility)
(Wx, Wy) - hydrodynamic body forces
Constitutive law and free boundary conditions:
actuated by time-dependent
preferred curvature  (s  ct)
EI - bending stiffness
 - viscoelastic damping
Rod shape tends to its preferred shape
 s - (actual) curvature
 - preferred curvature
Curvature defines shape.
s  

Case of time-independent preferred curvature  and no body forces:
  0 (no viscoelastic damping)
shape oscillates around preferred shape indefinitely
  0 (viscoelastic damping positive)
shape approaches preferred 
shape
Approximation of hydrodynamic forces
(Following G.I. Taylor*, …to avoid doing Navier-Stokes…)
a
v| |
v
v
v

Decompose forces in normal

and tangential components:
Normal forces
(neglecting drag)
N  CN a v 2
W = Nn + Lt
N  CN a v 2  8a v 3 / 2
Normal force proportional
to the product of diameter

and square of velocity.

2
L  2.7 2a v1/
 v ||
 - density,  - viscosity

214, 158-183, 1952
*Proc. Roy. Soc. Lond. A
Discretization of the rod: a chain of rigid links
Use finite differences in space s ih :

This mathematical discretization
has a nice physical interpretation
in terms of the segmented spinal
cords of eels and lampreys.
Discretization: springs, dashpots and muscles
Consider moments exerted at joints:
Force acting on
joint i:
From discretized
moments:
Ý
GR,L (t)  f R,L (t)   R,L  
R,L

We compute discrete
stiffnesses and curvatures:
Discretization: muscle properties
In the continuum (small h) limit the stiffness and curvature are:
fR  fL
EI  ab , and  
b
3
The dependence on material properties and body geometry is revealed.

Stiffness and curvature are now defined in terms
of body geomety, elastic properties and activation.
To complete the model we need to know what the
muscle forces are.
Neural activation and swimming in lamprey
From Fish and Wildlife.
The central pattern generator (CPG) of lamprey is a series of ipsi- and contralaterally
coupled neural oscillators distributed along the spinal notocord. In “fictive swimming”
in vitro, contralateral motoneurons burst in antiphase and there is a phase lag along
the cord from head to tail corresponding to about one full wavelength, at the typical 12 Hz burst frequency. This has been modeled as a chain of Kuramoto type coupled
rotators. The model can be justified by phase response and averaging theory:
[Cohen et al. J. Math Biol. 13, 345-369, 1982]
Incorporating muscle forces
Ý
GR,L (t)  f R,L (t)   R,L  
R,L

need to know this part
CD ≈ h - w i
d(CD)/dt ≈ -w di/dt
At each joint model the force on either side by muscle forces.
fR,L depends on:
(1) activation (calcium release, etc.) - traveling wave for now (CPG model?)
(2) length of muscle: h ± w i
(3) speed of muscle extension/contraction: ± w di/dt
A model of force development in lamprey
muscle
Williams, Bowtell and Curtin (*) developed a model for muscle forces based on a
simple kinetic model, using data obtained from isometric and ramp experiments.
The goal of this study was to construct a model of muscle tension development
which can reasonably predict the time course of muscle tension developed when
muscle is stimulated at different phases during sinusoidal movement, as occurs
during swimming.
The motivation for this study was to develop a model with adequate accuracy for
inclusion in a full neuromechanical model of the swimming lamprey.
(*) T.L. Williams, G. Bowtell, and N.A. Curtin. Predicting force
generation by lamprey muscle during applied sinusoidal movement
using a simple dynamic model. J. Exp. Biol. 201:869-875 (1998)
A. Peters & B. Mackay (1961). The structure
and innervation of the myotomes of the
lamprey.
J. Anat. 95, 575-585.
Muscle model chemical
constituents:
Output
of CPG
c: calcium ions
s: calcium-binding sites in
the sarcoplasmic
reticulum
f: calcium-binding sites in
the protein filaments
Mass action equations
d[c]/dt = k1[cs] - k2[c][s] - k3[c][f]
d[cf]/dt = k3[c][f] - k4[cf][f]
d[cs]/dt = -k1[cs] + k2[c][s]
d[f]/dt = -k3[c][f] + k4[cf][f]
d[s]/dt = k1[cs] - k2[c][s]
While the stimulus is on, k2=0.
While the stimulus is off, k1=0.
Reduced chemical kinetic equations
Constraints
[cs] + [c] + [cf] = CT
total # of calcium ions per litre is constant
[cs] + [s] = ST
total # of SR binding sites per litre is constant
[cf] + [f] = FT
total # of filament binding sites per litre is constant
5 equations in 5 variables plus 3 constraints
Variables: [c], [cf]
2 equations in 2 variables
Parameters: k1-k5, C, S, F
k1*(CT-[c]-[cf])
Stimulus on
k2[c](CT-ST-[c]-[cf])
Stimulus off
d[c]/dt = (k4*[cf]-k3[c])(FT-[cf]) +
d[cf]/dt = (k4*[cf]-k3[c])(FT-[cf])
Scaled Chemical Equations
All concentration variables and parameters are made nondimensionable by dividing by FT:
FT/FT = 1
Chosen ad hoc
CT/FT = C
C=2
Twice as much calcium is available
than needed to bind all the filaments.
ST/FT = S
[cf]/FT = Caf
thus Caf ≤ 1
[c]/FT = Ca
and Ca ≤ C
S=6
Thrice as many binding sites are
available in the SR than is required to bind all the
calcium.
k1*(C-Ca-Caf)
Stimulus on
dCa/dt = (k4*Caf-k3*Ca)(1-Caf) +
k2*Ca*(C-S-Ca-Caf) Stimulus off
dCaf/dt = (k4*Caf-k3*Ca)(1-Caf)
Mechanical model of muscle (A.V. Hill, 1938)
μS
LS
LC
L
L = L C + LS
TC = PC (Caf, LC, VC)
TS = μS * (LS - LS0) = PC
T P = μP * L
T = P C + TP
PC = T - TP
μP
LC(t) = L(t) - LS0 - PC(t)/μS
VC(t) = V(t) - (dPC/dt)/μS
Muscle properties
length-tension: force generated depends on muscle
length.
Investigate using isometric experiments.
force-velocity: force generated depends on speed of
lengthening or shortening of muscle.
Investigate using ramp experiments.
level of activation: force generated depends on the
number of muscle fibers activated and the frequency of that
activation.
Investigate by electrically stimulating muscle directly.
Basic assumptions of muscle model
1. The force developed is proportional to the number of calcium-activated
filaments.
2. Both the length-dependence and the velocity-dependence can be
described by independent multiplicative factors.
Pc = Pmax * Caf * λ(Lc) *
α (Vc)
Muscle experiments
lis and a are for a particular
lis
preparation
stimulating
electrode
a
output: force
required
Servo
motor
a. without
stimulation
b. with stimulation
l, dl/dt
input: desired length,
velocity
measure: length, velocity
Isometric experiments: constant muscle length
Ramp experiments: constant dl/dt
Sinusoidal experiments: l =
lis sin (wt)
Total measured force - passive force (mN)
Isometric tetanic contractions -- length dependence
P= Pmax * λ(Lc)
λ(Lc ) = 1 + λ2(Lc-Lc0)2
Total measured force - passive force (mN)
Ramp experiments -- velocity dependence
PC = Pmax * Caf * λ(LC) * α (VC)
αm * vc
vc < 0
αp * vc
vc ≥ 0
α (vc) = 1 +
Model equations
k1*(C-Ca-Caf)
Stimulus on
dCa/dt = (k4*Caf-k3*Ca)(1-Caf) +
k2*Ca*(C-S-Ca-Caf) Stimulus off
dCaf/dt = (k4*Caf-k3*Ca)(1-Caf)
Lc(t) = L(t) - LS0 - PC(t)/μS
Vc(t) = V(t) - (dPC/dt)/μS
PC = Pmax * Caf * λ(LC) * α (VC)
dP/dt = k5 * (PC - P)
necessary for fitting to data
Model parameters
Chosen ad hoc
C=2
Twice as much calcium is available than needed to bind all the
filaments.
S=6
Thrice as many binding sites are available in the SR than is
required to
bind all the calcium.
k5=100
Chosen large enough that P closely follows Pc.
Determined from the isotonic and ramp experiments:
αm
αp
λ2
Pmax
Found by least-squares fit to middle-length isometric data:
k1, k2, k3, k4
Sinusoidal experiments & predictions
Moment Dependence
Now we have that the moment depends on:
•Activation
•Curvature
•Rate of change of curvature
Depends on activation and state of the rod.
Swimming
Equal activations on both sides produces “straight” line swimming
QuickTime™ and a
decompressor
are needed to see this picture.
Swimming: Turns
Unequal activations on the sides produces turns.
QuickTime™ and a
decompressor
are needed to see this picture.
Shapes in time
Phase lags
lamprey
simulation
It’s qualitatively correct.
Comparison of effects
What’s happening
Summary
• Muscle model connected to rod “works”: it swims!
• Captures qualitatively the correct behavior (phase lags, shapes, etc.)
• Model allows flexibility to explore various effects
Future work
•More realistic fluid dynamics model (Navier-Stokes) & fluid-rod interaction
(Immersed Boundary Method)
•Better muscle model: Need effects of feedback and memory to get correct
isometric and dynamic behavior. Connect to proprioceptive and exteroceptive
sensing.
•Connect models of CPG, motoneurons, muscle force interaction, fluid dynamics
. . . “neurons to movement”
~~~~~~~~~~~~~~~~
References
•
•
•
•
•
•
•
•
G. Bowtell & T. Williams. Anguilliform body dynamics: Modeling the interaction
between muscle activation and body curvature. Phil. Trans. Roy. Soc. B 334:385-390 (1991)
A.H. Cohen, P. Holmes and R.H. Rand. The nature of coupling between segmental
oscillators of the lamprey spinal generator for locomotion. J. Math Biol. 13:345-369 (1982)
O. Ekeberg. A combined neuronal and mechanical model of fish swimming. Biol. Cyb.
69:363-374 (1992)
T. McMillen and P. Holmes. An elastic rod model for anguilliform swimming. J. Math.
Biol. 53:843-886 (2006)
T. McMillen, T. Williams and P. Holmes. Nonlinear muscles, viscoelastic damping and
body taper conspire to create curvature waves in the lamprey. In review, PLOS Comp. Biol.
G.I. Taylor. Analysis of the swimming of long and narrow animals. Proc. Roy. Soc. Lond. A
214:158-183 (1952)
T.L. Williams, G. Bowtell, and N.A. Curtin. Predicting force generation by lamprey
muscle during applied sinusoidal movement using a simple dynamic model. J. Exp. Biol.
201:869-875 (1998)
T.L. Williams, S. Grillner, V.V. Smoljaninov, P. Wallen and S. Rossignol. Locomotion in
lamprey and trout: The relative timing of activation and movement. J. Exp. Biol. 143:559566 (1989)