Predicting the effects of spatiotemporal modifications of muscle activation on the tentacle extension in squid

Squid use eight arms and two slender tentacles to capture prey. The muscular stalks of the tentacles are elongated approximately 80% in 20–40 ms towards the prey, which is adhered to the terminal clubs by arrays of suckers. Using a previously developed forward dynamics model of the extension of the tentacles of the squid Doryteuthis pealeii (formerly Loligo pealeii), we predict how spatial muscle-activation patterns result in a distribution of muscular power, muscle work, and kinetic and elastic energy along the tentacle. The simulated peak extension speed of the tentacles is remarkably insensitive to delays of activation along the stalk, as well as to random variations in the activation onset. A delay along the tentacle of 50% of the extension time has only a small effect on the peak extension velocity of the tentacle compared with a zero-delay pattern. A slight delay of the distal portion relative to the proximal has a small positive effect on peak extension velocity, whereas negative delays (delay reversed along stalk) always reduce extension performance. In addition, tentacular extension is relatively insensitive to superimposed random variations in the prescribed delays along the stalk. This holds in particular for small positive delays that are similar to delays predicted from measured axonal diameters of motor neurons. This robustness against variation in the activation distribution reduces the accuracy requirements of the neuronal control and is likely due to the non-linear mechanical properties of the muscular tissue in the tentacle.


SUMMARY OF THE MODEL
An overview of the definition of the symbols used for the model is provided in Supplementary Material: Table S1.Following van Leeuwen and Kier (1997), we assumed that the nominal tensile stress in the extensor muscles σ r (i.e., the tensile force per cross-sectional area of the initial relaxed state; subscript r refers to the radial direction) depends on the maximum active isometric stress σ max at optimum sarcomere length ℓ 0sarc , the normalized active state f a , the velocity dependence function f v , the filamentary overlap function f ℓ , and a passive component σ pas .
The subscript i refers to the segmental number.
The activation of each muscular element was described with a function f a,i in the interval [0, 1].The upper limit of the active state is described by f a,i = 1, and f a,i = 0 represents an inactive muscle.The activation is assumed to rise from zero to the maximum level, and then to remain maximal.Delays in the onset of the activation among tentacular segments can be set in an arbitrary manner.The activation function was described by the following equations: ) where t is time, t d,i is the activation delay with respect to the proximal segment, t a is the time between the start and full activation, and parameter q allows us to modify the basic sinusoidal shape.The activation along the tentacle is an essential input of the model, and the effects of its variation on the extension performance are the focus of this paper.
Following van Leeuwen and Kier (1997), we assumed that where ℓ myo,i and ℓ bz are the lengths of the myosin filament and the central bare zone on the myosin filament (which lacks myosin heads and cannot form cross-bridges).Furthermore, σ max,ref = 280 kPa and ℓ myo,ref = 1.58 µm represent values for a reference sarcomere.We assumed a fixed value of 0.14 µm for ℓ bz .
The non-linear functions f v , which depends on the longitudinal strain rate of the muscle fiber, and f ℓ , which depends on the strain (and hence the overlap in of the myosin and actin filaments in the sarcomere), are key properties of the muscle material.
The dependence of the active force on the non-dimensional strain rate ε * r,i was described by: where k is a constant and ε * r,i = εr,i / εmin,i ; εr,i is the strain rate and εmin,i is the minimum (unloaded) strain rate.Equation (S6), formulated by Otten (1987), is based on stretch experiments of vertebrate muscles Aubert (1956).Equation (S7) represents the equation of Hill (1938).
The filamentary overlap function was defined after van Leeuwen (1991): where, ℓ act,i is the summed length of two opposing actin filaments in one sarcomere, ℓ z is the width of the Z-disc, ℓ sarc,i is the sarcomere length, and D act and D myo account for cross-bridge losses due to actin overlap and interaction between myosin filament and Z-disc.Finally, C myo accounts for resistive forces owing to the collision of the myosin filaments with the Z-disc of the sarcomere.The reference sarcomere length for active force production is defined as: All sarcomeres were prescribed to be at ℓ 0sarc in the initial relaxed state of the tentacle and in each segment all sarcomeres were assumed to have the same length.Hence, ℓ sarc,i was calculated from: The passive component in equation ( S1) was expressed as (van Leeuwen and Kier, 1997): where ε c is the critical strain above which the relationship is linear, and c 1pas , • • • , c 4pas are constants.ε i stands either for ε ℓ,i or for ε r,i .To ensure a continuous first derivative at ε c , we prescribed c 1pas and c 2pas as: Following van Leeuwen and Kier (1997), the passive stiffness of the longitudinal muscle fiber bundles was assumed to be equal to that of the extensor muscles.
From the requirement of constant volume in each tentacle segment, van Leeuwen and Kier (1997) and for the relationships between the radial velocity ṙi and radial acceleration ri of segmental peripheral boundary i, and its radius r i , its height h i , and the time derivatives of its height ḣi and ḧi .
van Leeuwen and Kier (1997) derived for the pressure p i in tentacle segment i where η represent a fraction (taken 0.7 at the base and in most cases 0.6 at distal end of the stalk to allow for muscular tapering) to account for the volume occupied by the extensor muscles in the segment, A s,i is the instantaneous outer surface of the segment, A s0,i is the outer surface at rest, m i is the mass of segment i, and α is a constant for defining the added mass of water.Values of η were linearly interpolated between those at base and tip of the stalk.At the tip, we also used η = 0.4, 0.5 and 0.7 to examine the effect of muscle tapering on extension performance.
The force exerted by segment i on boundary i (located at the distal side of the segment), was computed by van Leeuwen and Kier (1997) as: where A c,i is the cross-sectional area of disc i, µ is a viscous friction parameter associated with the strain rate εℓ,i in the longitudinal muscle fiber bundles, and A cℓ0 and σ ℓ,i are, respectively, the reference cross-sectional area and longitudinal stress of the longitudinal muscle fiber bundles.The external force at the most distal boundary of the club is prescribed as zero: Following van Leeuwen and Kier (1997), we modelled the tentacular club as a single rigid segment.
We can write the acceleration a i of boundary b i as: where m b,i is the mass at the boundary between segments i and i + 1.Expression (S24) is an implicit equation for a i , because F i and F i+1 depend on ri and ri+1 , whereas, in equation (S20), ḧi = a i − a i−1 .However, van Leeuwen and Kier (1997) showed that at each instant, the accelerations of the segment boundaries a 1 • • • a n−1 can be obtained by solving the following matrix equation: Expression (S25) is a tridiagonal system of linear equations because each boundary mass has direct mechanical interactions only with its proximal and distal neighbors (see force diagram of Figure 2c in the main text).The components c change with time.Expressions for the components c in equation (S25) can derived as follows.First, we define: By substitution, equations (S20), (S21), and (S22) can now be written as: By subsequent substitution of (S32) into (S33), ( S33) into (S34), and (S34) into (S24), the following equation is obtained: where The c components in equation ( S25) can now be derived from these equations by substituting appropriate values for i, while taking into account the boundary conditions at the base and the tip of the stalk.For the base, we prescribed a constant forward velocity.For the tip, the pressure was set to zero.A hydrodynamic pressure acts at the tentacular club, which was simulated with an added mass term (see van Leeuwen and Kier (1997)).

Effect of deterministic activation delays and additive noise on peak-extension velocity
Figure 3b (main document) shows the probability density distribution that we used for the simulations of the noise modulated motor input delays along the tentacular stalk.Table S2 shows (1) the peak extension velocities (v max ) of the stalk for the deterministic delay cases −10, −5, 0, 5, 10, and 15 ms and (2) the median, mean and standard deviations of Monte Carlo simulations (n = 100 for each case) of the effects of the applied additive noise in the timing of the motor input on the maximum extension velocity for the probability intervals of [−5 ms 5 ms] and [−10 ms 10 ms], with a single additive random delay per segment.Table S3 shows the effects for the same deterministic delay cases and noise intervals, but now ten additive noise samples were used per segment to compute the weighed activation for the overall segmental activation.
Table S2.Median, mean, and standard deviation (Std) of the maximum extension velocity, v max , of the tentacular stalk for a simulated series of six different deterministic delays, ∆t d , with an added noise distribution in the range [−5 ms 5 ms] and [−10 ms 10 ms] (corresponding respectively to Figures 11 and 12 in the main text) of the motor input along the stalk.Per segment, a single random delay number was picked from the beta distributed probability density function (see Figure 3b).One-hundred runs were made for each noise-modulated deterministic delay (n = 100).The first data row shows the maximum extension velocities without motor-input noise.S3.Median, mean, and standard deviation (Std) of the maximum extension velocity, v max , of the tentacular stalk for a simulated series of six different deterministic delays, ∆t d , with an added noise distribution in the range [−5 ms 5 ms] and [−10 ms 10 ms] (corresponding respectively to Figures 13 and 14 in the main document) of the motor input along the stalk.Per segment, ten random delays were picked from the beta-distributed probability density function (see Figure 3b).The segmental random delays were averaged so as to generate a single active state function for each segment.One-hundred runs were made for each noise-modulated deterministic delay (n = 100).The first data row shows the maximum extension velocities without motor-input noise.

Bayesian statistics
We used the Bayesian estimation approach of Kruschke (2013) to compare (1) the mean peak extension performance of additive noise affected tentacle extensions with the corresponding v max for each of the six deterministic cases (i.e., −10, −5, 0, 5, 10, and 15 ms) and (2) the mean peak extension performances of the simulated per-segment additive noise affected deterministic activation delays along the tentacular stalk (for two different noise windows, single or average segmental noise perturbations, and six different values of ∆t d , this leads to a total of 4 times 15 = 60 comparisons).We adapted the MATLAB toolbox for Bayesian estimation of Winter (2016) to accommodate an automatized analysis of an array of two-group comparisons.We verified the validity of the code by analyzing the data set used by Kruschke (2013) to demonstrate the approach.We used the following settings of the input parameters for the Bayesian analysis: number of saved steps = 50000; number of chains = 4; number of thinning steps = 1; number of burn-in steps = 2000.These settings guarantee a robust Bayesan estimation (see Kruschke (2013) for a discussion of the settings of these parameters).
Figure S1 illustrates a Bayesian estimation example of the comparison between the additive noise affected activation delay cases −10 ms and 15 ms, with noise probability interval [−5 ms 5 ms]).Interestingly, the analysis shows that for these extreme negative and positive delay cases the mean peak extension velocities are greater than the corresponding deterministic v max values.
Similarly, Figure S2 illustrates a Bayesian estimation example with the comparison between the noise affected activation delay cases 0 ms and 5 ms, with noise probability interval [−10 ms 10 ms].For these cases, the mean extension velocities show a substantial decrease compared with the corresponding deterministic v max (order 10 % reduction).The noise disturbed ∆t d = 5 ms yields a significantly higher mean v max than the equivalent ∆t d = 0 ms case, which is in line with the deterministic performance difference.

Figure S2 .
Figure S2.Distributions of credible values generated by the Bayesian estimation procedure.Number of input simulation samples per group is 100.Labeling is as in Figure S1.(a) Distribution for the mean of v max for the per-segment noise disturbed zero activation delay along the tentacular stalk.The green vertical line shows the deterministic v max = 2.254 m s −1 .The dashed vertical lines indicate the boundaries of the velocity ROPE, [v max − 0.01 m s −1 v max + 0.01 m s −1 ].All credible values are outside the ROPE and smaller than the deterministic v max , demonstrating a negative effect of the additive noise on the peak extension velocity.(b) Idem for the corresponding standard deviation.(c) Distribution for the mean of v max for the per-segment noise disturbed activation delay of 5 ms along the tentacular stalk.(d) Idem for the associated standard deviation.(e) Distribution for the difference of the means of the two groups.All credible values are below zero and located outside the ROPE.This indicates that the mean peak extension velocity of the noise perturbed ∆t d = 5 ms case is greater than that of the equivalent perturbed zero delay case.( f ) Idem for the difference of the standard deviations.(g) Distribution for the log-transformed shape parameter ν.(g) Distribution for the effect size with none of the values in the ROPE, [−0.01 0.01].This underpins the conclusion that the perturbed ∆t d = 5 ms case performs better than the equivalent 0 ms case.

Table S1 :
Symbols, definitions, and assigned values  (continued)εmin,ref = −17 s −1 , which is equivalent to V max = 17 ℓ s −1 η ivolume fraction of extensor muscles in segment i. η is linearly interpolated between the proximal and distant segments of the tentacular stalk; η 1 = 0.7 for all simulations ∆t d deterministic activation delay between the tip and base segments of the tentacular stalk.∆t d,vmax the ∆t d that results in the highest peak extension velocity of the tentacular stalk.
8 gram m i mass of tentacle segment i; m i = πr 2 0 h 0,i ρ n number of segments of tentacle model; n = 51 (50 stalk segments, 1 club segment); n can also refer to the number of data in a group x variable in t distribution (equation (2) in the main text) α constant for defining the added mass of water; α = 0.08