In collaboration with Thelma Williams and Philip Holmes.

1. Nonlinear muscles, viscoelasticity and body taper in the creation of curvature wavesSIAM PDEs December 10, 2007~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Tyler McMillenCalifornia State University, Fullerton In collaboration with Thelma Williams and Philip Holmes.

2. Relative timing of activation and movement Curvature travels slower than activation. Figure from Williams, et. al., J. Exp. Biol. (1989)

3. How are waves of curvature created and propagated? Why does curvature travel slower than activation?

4. Outline • Elastic rod model • Resistive fluid forces • Discretization -- chain of rigid links • Muscle forces

5. g Dynamics of an actuated elastic rod in the plane 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 EI - bending stiffness - viscoelastic damping

6. Rod shape tends to its preferred shape Curvature defines shape. Case of time-independent preferred curvature and no body forces: (no viscoelastic damping) shape oscillates around preferred shape indefinitely (viscoelastic damping positive) shape approaches preferred shape

7. Approximation of hydrodynamic forces(Following G.I. Taylor*, …to avoid doing Navier-Stokes…) a v v Decompose forces in normal and tangential components: Normal forces (neglecting drag) W = Nn + Lt Normal force proportional to the product of diameter and square of velocity. - density, - viscosity *Proc. Roy. Soc. Lond. A 214, 158-183, 1952

8. Discretization of the rod: a chain of rigid links Use finite differences in space : This mathematical discretization has a nice physical interpretation in terms of the segmented spinal cords of eels and lampreys.

9. Discretization: springs, dashpots and muscles Consider moments exerted at joints: Force acting on joint i: From discretized moments: We compute discrete stiffnesses and curvatures:

10. Discretization: muscle properties In the continuum (small h) limit the stiffness and curvature are: 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.

11. 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 1-2 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]

12. Incorporating muscle forces 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

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

14. A. Peters & B. Mackay (1961). The structure and innervation of the myotomes of the lamprey. J. Anat. 95, 575-585.

16. Output of CPG Muscle model chemical constituents: c: calcium ions s: calcium-binding sites in the sarcoplasmic reticulum f: calcium-binding sites in theprotein filaments

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

18. Reduced chemical kinetic equations 5 equations in 5 variables plus 3 constraints 2 equations in 2 variables Variables: [c], [cf] Parameters: k1-k5, C, S, F 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 k1*(CT-[c]-[cf]) Stimulus on d[c]/dt = (k4*[cf]-k3[c])(FT-[cf]) + k2[c](CT-ST-[c]-[cf]) Stimulus off d[cf]/dt = (k4*[cf]-k3[c])(FT-[cf])

19. All concentration variables and parameters are made non-dimensionable by dividing by FT: FT/FT = 1 CT/FT = C ST/FT = S [cf]/FT = Caf thus Caf ≤ 1 [c]/FT = Ca and Ca ≤ C Scaled Chemical Equations 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. 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)

20. L Mechanical model of muscle (A.V. Hill, 1938) μS LS LC μP L = LC + LS TC = PC (Caf, LC, VC) TS= μS* (LS - LS0) = PC TP = μP*L T = PC + TP PC= T - TP LC(t) = L(t)- LS0 - PC(t)/μS VC(t) = V(t) - (dPC/dt)/μS

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

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

23. Muscle experiments output: force required a. without stimulation b. with stimulation stimulating electrode Servo motor l, dl/dt measure: length, velocity Isometric experiments: constant muscle length Ramp experiments: constant dl/dt Sinusoidal experiments: l = lis sin (wt) lis and a are for a particular preparation lis a input: desired length, velocity

24. Isometric tetanic contractions -- length dependence P= Pmax * λ(Lc) Total measured force - passive force (mN) λ(Lc ) = 1 + λ2(Lc-Lc0)2

25. Ramp experiments -- velocity dependence Total measured force - passive force (mN) PC = Pmax * Caf * λ(LC) * α (VC) αm * vc vc < 0 α (vc) = 1 + αp * vc vc ≥ 0

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

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

28. Sinusoidal experiments & predictions

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

30. Swimming Equal activations on both sides produces “straight” line swimming

31. Swimming: Turns Unequal activations on the sides produces turns.

32. Shapes in time

33. Phase lags lamprey simulation It’s qualitatively correct.

34. Comparison of effects

35. What’s happening

36. 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” ~~~~~~~~~~~~~~~~

37. 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:559-566 (1989)