Control for multifunctionality: bioinspired control based on feeding in Aplysia californica

Abstract

Animals exhibit remarkable feats of behavioral flexibility and multifunctional control that remain challenging for robotic systems. The neural and morphological basis of multifunctionality in animals can provide a source of bioinspiration for robotic controllers. However, many existing approaches to modeling biological neural networks rely on computationally expensive models and tend to focus solely on the nervous system, often neglecting the biomechanics of the periphery. As a consequence, while these models are excellent tools for neuroscience, they fail to predict functional behavior in real time, which is a critical capability for robotic control. To meet the need for real-time multifunctional control, we have developed a hybrid Boolean model framework capable of modeling neural bursting activity and simple biomechanics at speeds faster than real time. Using this approach, we present a multifunctional model of Aplysia californica feeding that qualitatively reproduces three key feeding behaviors (biting, swallowing, and rejection), demonstrates behavioral switching in response to external sensory cues, and incorporates both known neural connectivity and a simple bioinspired mechanical model of the feeding apparatus. We demonstrate that the model can be used for formulating testable hypotheses and discuss the implications of this approach for robotic control and neuroscience.

Appendices

Appendices

1.1 Semi-implicit integration scheme

Suppose a continuously varying quantity x satisfies the initial value problem

$$\begin{aligned} \frac{\mathrm{d}x(t)}{\mathrm{d}t}=\frac{-(x(t)-x_\infty (y_2(t),\ldots ,y_n(t)))}{\tau },\quad x(t_0)=x_0 \end{aligned}$$

(11)

where $x_\infty (t)$ is set by the other variables in our system, say $\{y_i\}_{i=2}^n$, generally following some nonlinear dependencies, and $\tau$ is a fixed time constant. We would like to implement a numerical approximation to the exact solution for x, namely

$$\begin{aligned} x(t)=e^{-(t-t_0)/\tau }x_0+\frac{1}{\tau }\int _{t_0}^t e^{-(t-s)/\tau }x_\infty (y_2(s),\ldots ,y_n(s))\,\mathrm{d}s, \end{aligned}$$

(12)

along with the remaining variables that satisfy their own differential equations. Euler’s forward method is convenient to implement but prone to numerical instability. Euler’s backward or implicit method is numerically stable but computationally expensive, as it requires solving an implicit equation at each step. Both methods proceed from a discrete approximation of the derivative, namely

$$\begin{aligned} \frac{x(t+h)-x(t)}{h}\approx \frac{\mathrm{d}x}{\mathrm{d}t}. \end{aligned}$$

(13)

In both cases we create an update rule $x(t)\rightarrow x(t+h)$, by evaluating the right-hand side of (13) at either time t or time $t+h$, and solving for $x(t+h)$.

Forward:

$$\begin{aligned} \frac{x(t+h)-x(t)}{h}&=\frac{ -(x(t)-x_\infty (y_2(t),\ldots ,y_n(t))}{\tau }\end{aligned}$$

(14)

$$\begin{aligned} x(t+h)&=x(t)-\frac{x(t)-x_\infty (y_2(t),\ldots ,y_n(t)}{\tau }h \end{aligned}$$

(15)

Backward:

$$\begin{aligned}&\frac{x(t+h)-x(t)}{h}\nonumber \\&\quad = \frac{-(x(t+h)-x_\infty (y_2(t+h),\ldots ,y_n(t+h))}{\tau } \end{aligned}$$

(16)

$$\begin{aligned} x(t+h)&=\frac{\tau x(t)+h x_\infty (y_2(t+h),\ldots ,y_n(t+h))}{\tau +h}. \end{aligned}$$

(17)

Since the variables $y_2,\ldots ,y_n$ appear on the right-hand side of (17) evaluated at the later time point, $t+h$, (17) is part of a system of n nonlinear equations that must be solved simultaneously to determine the system state at $t+h$. Both numerical schemes (15) and (17) are first-order accurate, meaning that the truncation error between the true solution (12) and the numerical approximation scales as $\mathcal {O}(h^2)$ on each time step, with a global error (after T/h time steps for a simulation of total runtime T) that is $\mathcal {O}(h)$.

Semi-implicit: In our model implementation, we use a semi-implicit method based on the approximation

$$\begin{aligned} \frac{x(t+h)-x(t)}{h}\approx \frac{-(x(t+h)-x_\infty (y_2(t),\ldots ,y_n(t))}{\tau }, \end{aligned}$$

(18)

namely

At each time step we update x using a weighted average of its past value x(t) and its target value $x_\infty (t)$, with the (short) timestep h and the intrinsic time constant $\tau$ providing the relative weight of past and future. We expect an accurate approximation to (12) provided $h\ll \tau$. As we show below, the method is first-order accurate, and numerically stable, but it does not require solving an implicit equation at each time step. Thus, this method combines the advantages of both the forward and backward methods. The method may be seen as an example of operator splitting.[139]

To see that (19) is first-order accurate, we assume that x(t) is smooth enough to have Taylor expansions through the second order. Thus, for $h\ll 1$ we may write

$$\begin{aligned} x(t+h)&=x(t)+h\frac{\mathrm{d}x}{\mathrm{d}t}(t)+\mathcal {O}(h^2),\quad \text { as }h\rightarrow 0 \end{aligned}$$

(20)

$$\begin{aligned}&=x(t)+h\frac{x_\infty (\mathbf {y}(t))-x(t)}{\tau }+\mathcal {O}(h^2) \end{aligned}$$

(21)

$$\begin{aligned}&=x(t)\frac{\tau +h}{\tau +h}+h\frac{x_\infty (\mathbf {y}(t))-x(t)}{\tau +h}\left( \frac{\tau +h}{\tau }\right) +\mathcal {O}(h^2) \end{aligned}$$

(22)

$$\begin{aligned}&=\frac{\tau x(t)}{\tau +h}+\frac{hx(t)}{\tau +h}+\left( \frac{\tau +h}{\tau }\right) \frac{hx_\infty (\mathbf {y}(t))}{\tau +h}\nonumber \\&\quad -\left( \frac{\tau +h}{\tau }\right) \frac{hx(t)}{\tau +h}+\mathcal {O}(h^2) \end{aligned}$$

(23)

$$\begin{aligned}&=\frac{\tau x(t)}{\tau +h}+\frac{hx_\infty (\mathbf {y}(t))}{\tau +h}+\frac{hx(t)}{\tau +h}-\frac{hx(t)}{\tau +h}+\mathcal {O}(h^2) \end{aligned}$$

(24)

$$\begin{aligned}&=\frac{\tau x(t)}{\tau +h}+\frac{hx_\infty (\mathbf {y}(t))}{\tau +h}+\mathcal {O}(h^2),\quad \text { as }h\rightarrow 0. \end{aligned}$$

(25)

Thus, the semi-implicit scheme (19) is first-order accurate in the time step h.

To see that (19) is numerically stable, suppose that we fix $\mathbf {y}$ so that $x_\infty (\mathbf {y})=c,$ a constant. Clearly if $x(t)=c$ then $x(t+h)=c$ as well, so $x=c$ is a fixed point of the iteration (19), under this assumption. Numerical stability follows if we can show that $x=c$ is a stable fixed point for all $h>0$, as we now establish. Let $x(t_0+nh)=c+a_n,$ with $a_0$ arbitrary. Then

$$\begin{aligned} a_{n+1}&=x(t_0+nh+h)-c \end{aligned}$$

(26)

$$\begin{aligned}&=\frac{\tau x(t_0+nh)+hc}{\tau +h}-c \end{aligned}$$

(27)

$$\begin{aligned}&=\frac{\tau (c+a_n)+hc}{\tau +h}-c \end{aligned}$$

(28)

$$\begin{aligned}&=\frac{\tau }{\tau +h}a_n\rightarrow 0,\quad \text { as }n\rightarrow \infty \end{aligned}$$

(29)

no matter the size of the timestep $h>0$. Thus, the scheme (19) is both (first-order) accurate and numerically stable.

The head and grasper position variables $x_\text {h}$, $x_\text {g}$ form a linearly coupled pair, for which we can extend the semi-implicit algorithm given in one-dimensional form above. In general, consider a nonhomogeneous linear system expressed in terms of a vector $\mathbf {x}$, a matrix A, and a forcing vector $\mathbf {b}$:

$$\begin{aligned} \frac{\mathrm{d}\mathbf {x}}{\mathrm{d}t}=A(t)\mathbf {x}(t)+\mathbf {b}(t). \end{aligned}$$

(30)

To set up a semi-implicit first-order iteration scheme, observe that

$$\begin{aligned}&\frac{\mathbf {x}(t+h)-\mathbf {x}(t)}{h}=A(t)\mathbf {x}(t+h)+\mathbf {b}(t)+\mathcal {O}(h^2),\quad \text { so } \end{aligned}$$

(31)

$$\begin{aligned}&\mathbf {x}(t+h)-hA(t)\mathbf {x}(t+h)=\mathbf {x}(t)+h\mathbf {b}(t)+\mathcal {O}(h^2),\quad \text { therefore } \end{aligned}$$

(32)

$$\begin{aligned}&\mathbf {x}(t+h)=\left( I-hA(t)\right) ^{-1}(\mathbf {x}(t)+h\mathbf {b}(t))+\mathcal {O}(h^2). \end{aligned}$$

(33)

Dropping the $\mathcal {O}(h^2)$ term and writing the update scheme in MATLAB style notation gives

(Here the backslash notation $M\backslash \mathbf {u}$ stands for $M^{-1}\mathbf {u}$, i.e., the least-squares solution $\mathbf {y}$ to the linear system $M\mathbf {y}=\mathbf {u}$.) Comparing this update scheme with (19), it is easy to check that they are consistent for a single variable by setting $A=-1/\tau$ and $b=x_\infty /\tau$.

Table 3 Table of symbols

In our case the head and grasper positions $x_\text {h}$ and $x_\text {g}$ comprise a linearly coupled system, $\mathbf {x}$, and the coupling matrix A is $2\times 2$, so $I-hA$ can be inverted explicitly, provided h is smaller than the reciprocal of the largest positive eigenvalue of A. (If no eigenvalues of A are positive real numbers, then $I-hA$ can always be inverted.) For a general $2\times 2$ system the update rule reads

$$\begin{aligned}&\left( \begin{array}{cc}x_1(t+h) \\ x_2(t+h)\end{array}\right) = \frac{1}{1-h\text {Tr}A +h^2\det A}\nonumber \\&\quad \quad \left( \begin{array}{c} (1-hA_{22})(x_1+h b_1)+hA_{12}(x_2+hb_2)\\ hA_{21}(x_1+hb_1)+(1-hA_{11})(x_2+hb_2) \end{array} \right) , \end{aligned}$$

(35)

with all time-varying elements of the right-hand side evaluated at time t. In (35) $\text {Tr}A$ and $\det A$ denote the trace and determinant of A, respectively. Thus, truncating terms of order $\mathcal {O}(h^2)$ and higher gives a first-order semi-implicit update scheme for two-component state and forcing vectors $\mathbf {x}$ and $\mathbf {b}$:

In (36) I denotes the $2\times 2$ identity matrix.

1.2 Table of symbols

Table 3 provides a table of symbols.

1.3 Boolean logic of Aplysia feeding control

The following sections detail the logic implementations for the activity of each neuron in the controller network. All neurons, except for B4/B5, are implemented as standard Boolean elements and are either OFF (0), or ON (1). B4/B5 is a ternary unit which is either OFF (0), weakly ON (1), or strongly ON (2). The activity state of B4/B5 must be checked using conditional logic prior to negation, as $\,!\,N_{B4/B5}$ is undefined (see Sect. 3.5). These interactions include known direct connections between neurons based on previous literature as well as hypothesized connections that may be direct or indirect, and indirect sensory feedback pathways. Sensory feedback pathways involving proprioception of the grasper position and pressure exerted by the closed grasper are gated by logic tests of position and pressure relative to user-specified thresholds. Allowing these thresholds to vary in response to activity levels of interneurons and external sensory cues provides an approximation of neuromodulation. All time varying elements on the right side of the logic equations are at time (j).

Cerebral Interneurons

• 1. Metacerebral Cell

  $$\begin{aligned} N_\text {MCC}(j+1) = [\text {arousal}] \end{aligned}$$

  (37)

• 2. CBI-2

  CBI-2 is activated by sensory inputs present in biting and rejection, but not in swallowing.

  $$\begin{aligned} \begin{aligned}&N_\text {CBI-2}(j+1) \\&\quad = N_\text {MCC}\,\,(\,!\,N_\text {B64}) \,\,\big (\left( [\text {lips}_\text {mech}]\,\,[\text {lips}_\text {chem}]\,\,\,!\,[\text {grasper}_\text {mech}]\right) \parallel \\&\qquad \left( [\text {grasper}_\text {mech}]\,\,\,!\,[\text {lips}_\text {chem}]\right) \big ) \end{aligned} \end{aligned}$$

  (38)

  With the hypothesized connections in Sect. 4.3, these equations change to:

  $$\begin{aligned} \begin{aligned}&N_\text {CBI-2}(j+1) \\&\quad = N_\text {MCC}\,\,(\,!\,N_\text {B64}) \,\,\big (\left( [\text {lips}_\text {mech}]\,\,[\text {lips}_\text {chem}]\,\,\,!\,[\text {grasper}_\text {mech}]\right) \parallel \\&\qquad \left( [\text {grasper}_\text {mech}]\,\,\,!\,[\text {lips}_\text {chem}]\right) \parallel \left( N_\text {B4/B5}\ge 2 \right) \big ) \end{aligned} \end{aligned}$$

  (39)

• 3. CBI-3

  CBI-3 is activated by sensory inputs present in biting and swallowing, but not in rejection.

  $$\begin{aligned} \begin{aligned} N_\text {CBI-3}(j+1) =&N_\text {MCC}\,\,[\text {lips}_\text {mech}]\,\,[\text {lips}_\text {chem}]\end{aligned} \end{aligned}$$

  (40)

  With the equations and refractory period proposed in Sect. 4.3, the logic implementation for CBI-3 changes to include a gating state variable based on whether or not the neuron is in a refractory state following strong inhibition. Similar logic could be added to other nodes in the network as needed based on animal experiments. This period was included here as part of the hypothesis that strong activation of B4/B5 triggers rejection in animals that are swallowing. This hypothesis and an assessment of whether this refractory period occurs in CBI-3 in animal preparations or whether this effect is due to another mechanism could be tested experimentally. The equation becomes:

  $$\begin{aligned} \begin{aligned} N_\text {CBI-3}(j+1) = \,&N_\text {MCC}\,\,[\text {lips}_\text {mech}]\,\,[\text {lips}_\text {chem}]\\&\,\,\left( N_\text {B4/B5}< 2 \right) \,\,\left( \,!\,[\text {refractory}_\text {CBI-3}]\right) \end{aligned} \end{aligned}$$

  (41)

• 4. CBI-4

  CBI-4 is activated by sensory inputs present in swallowing and rejection, but not in biting.

  $$\begin{aligned} N_\text {CBI-4}(j+1) = N_\text {MCC}\,\,\left( [\text {lips}_\text {mech}]\parallel [\text {lips}_\text {chem}]\right) \,\,[\text {grasper}_\text {mech}]\end{aligned}$$

  (42)

Buccal Interneurons

• 1. B64

  Activity in $N_\text {B64}$ is influenced by the activity of the $N_\text {MCC}$ and $N_\text {B31/B32}$. It is also excited by protraction and inhibited by retraction. The proprioceptive feedback is implemented as:

  $$\begin{aligned}&\text {B64}_\text {proprioception} \nonumber \\&\quad = (N_\text {CBI-3}\,\,(([\text {grasper}_\text {mech}]\,\,[\text {protracted}_{N_\text {B64},\text {swallow}}]) \parallel \nonumber \\&\qquad ((\,!\,[\text {grasper}_\text {mech}])\,\,[\text {protracted}_{N_\text {B64},\text {bite}}]))\quad ) \parallel \nonumber \\&\qquad ((\,!\,N_\text {CBI-3})\,\,[\text {protracted}_{N_\text {B64},\text {reject}}]) \end{aligned}$$

  (43)

  where

  $$\begin{aligned}&\left[ \text {protracted}_{N_\text {B64},\text {swallow}}\right] = [x_{\text {g/h}}> z_{\text {B64},\text {swallow}}] \end{aligned}$$

  (44)
  $$\begin{aligned}&\left[ \text {protracted}_{N_\text {B64},\text {bite}}\right] = [x_{\text {g/h}}> z_{\text {B64},\text {bite}}] \end{aligned}$$

  (45)
  $$\begin{aligned}&\left[ \text {protracted}_{N_\text {B64},\text {reject}}\right] = [x_{\text {g/h}}> z_{\text {B64},\text {reject}}] \end{aligned}$$

  (46)

  This amounts to the threshold being depend on the behavior with different threshold values for bites, swallows, and rejections.

  $$\begin{aligned} \begin{aligned} N_\text {B64}(j+1) =&N_\text {MCC}\,\,(\,!\,N_\text {B31/B32}) \,\,\text {B64}_\text {proprioception} \end{aligned} \end{aligned}$$

  (47)

• 2. B4/B5

  $N_\text {B4/B5}$ has been shown to have varying effects when firing strongly versus weakly. To represent this in the modeling framework, quiescence is represented as 0, weak firing as 1 and strong firing as 2. The neurons are quiescent during biting, and they fire weakly during the retraction phase of swallowing. The neurons fire strongly when stimulated with the external electrode and during the retraction phase of rejection. During rejection, B4/B5 is observed to cease firing, allowing B3/B6/B9 to fire briefly at the end of the behavior. To implement this, we have used a proprioceptive feedback pathway which inhibits the activity of B4/B5 once the grasper has reached a user-specified level of retraction.

  $$\begin{aligned} N_\text {B4/B5}(j+1)&= N_\text {MCC}\,\,\bigg ( (\,!\,[\text {electrode}_\text {B4/B5}]) \,\,\big ( 2 (\,!\,N_\text {CBI-3}) \,\,N_\text {B64}\nonumber \\&\quad \,\,[\text {protracted}_{N_\text {B4/B5}}] \nonumber \\&\quad +N_\text {CBI-3}\,\,[\text {grasper}_\text {mech}]\,\,N_\text {B64}\big ) \nonumber \\&\quad +2 \, [\text {electrode}_\text {B4/B5}] \bigg ) \end{aligned}$$

  (48)

  where

  $$\begin{aligned}{}[\text {protracted}_{N_\text {B4/B5}}] = [x_{\text {g/h}}> z_{\text {B4/B5}}] \end{aligned}$$

  (49)

• 3. B20

  $$\begin{aligned} \begin{aligned}&N_\text {B20}(j+1) \\&\quad = N_\text {MCC}\big (N_\text {CBI-2}\parallel N_\text {CBI-4}\parallel N_\text {B31/B32}\big ) \,\,\,!\,N_\text {CBI-3}\,\,\,!\,N_\text {B64}\end{aligned} \end{aligned}$$

  (50)

• 4. B40/B30

  $N_\text {B40/B30}$ has fast inhibitory and slow excitatory connections to $N_\text {B8}$. To capture this, we record the time (j) at which $N_\text {B40/B30}$ transitions between states for later use in the $N_\text {B8}$ activity calculations (see below). First, the activity of $N_\text {B40/B30}$ in the next time step is determined:

  $$\begin{aligned} \begin{aligned} N_\text {B40/B30}(j+1) = N_\text {MCC}\,\,\big (N_\text {CBI-2}\parallel N_\text {CBI-4}\parallel N_\text {B31/B32}\big ) \,\,\,!\,N_\text {B64}\end{aligned} \end{aligned}$$

  (51)

  After calculating the new activity, we assess transitions as defined by the following pseudocode:

  if $(N_\text {B40/B30}(j) == 0 \, \text {AND} \, N_\text {B40/B30}(j+1) == 1)$, then set $t_{N_\text {B40/B30}\text {,on}}$ = j;

  if $(N_\text {B40/B30}(j) == 1 \, \text {AND} \, N_\text {B40/B30}(j+1) == 0)$, then set $t_{N_\text {B40/B30}\text {,off}}$ = j;

Buccal Motor Neurons

• 1. B31/B32

  $N_\text {B31/B32}$ receives input from interneurons and proprioceptive feedback. To capture possible modulation of $N_\text {B31/B32}$ and generate multifunctional behavior under different sensory cues, behavior-dependent proprioceptive inputs are implemented. Though the resulting full equation for $N_\text {B31/B32}$ activity is large, it can be broken down to three sections: (1) if $N_\text {CBI-3}$ is active and there is sensory stimuli in the grasper (swallowing), (2) if $N_\text {CBI-3}$ is active and there is NOT sensory stimuli in the grasper (biting), and (3) if $N_\text {CBI-3}$ is NOT active (rejection).

  $$\begin{aligned}&N_\text {B31/B32}(j+1) \nonumber \\&\quad = N_\text {MCC}\,\,\bigg ( N_\text {CBI-3}\,\,\nonumber \\&[\text {grasper}_\text {mech}]( (\,!\,N_\text {B64}) \,\,((![\text {pressure}_{N_\text {B31/B32},\text {ingestion}}]) \parallel N_\text {CBI-2}) \,\,\nonumber \\&\quad ((\,!\,N_\text {B31/B32})\,\,[\text {retracted}_{N_\text {B31/B32},\text {swallow},\text {off}}] +\nonumber \\&\quad N_\text {B31/B32}\,\,[\text {retracted}_{N_\text {B31/B32},\text {swallow},\text {on}}]))\nonumber \\&(\,!\,[\text {grasper}_\text {mech}]) ( (\,!\,N_\text {B64}) \,\,((![\text {pressure}_{N_\text {B31/B32},\text {ingestion}}]) \parallel N_\text {CBI-2}) \,\,\nonumber \\&\quad ((\,!\,N_\text {B31/B32})\,\,[\text {retracted}_{N_\text {B31/B32},\text {bite},\text {off}}] +\nonumber \\&\quad N_\text {B31/B32}\,\,[\text {retracted}_{N_\text {B31/B32},\text {bite},\text {on}}])) + \nonumber \\&(\,!\,N_\text {CBI-3}) \,\,\big ((\,!\,N_\text {B64}) \,\,[\text {pressure}_{N_\text {B31/B32},\text {rejection}}] \,\,\nonumber \\&\quad (N_\text {CBI-2}\parallel N_\text {CBI-4}) \,\,\nonumber \\&\quad ((\,!\,N_\text {B31/B32})\,\,[\text {retracted}_{N_\text {B31/B32},\text {reject},\text {off}}] + \nonumber \\&\quad N_\text {B31/B32}\,\,[\text {retracted}_{N_\text {B31/B32},\text {reject},\text {on}}] \big ) \bigg ) \end{aligned}$$

  (52)

  where

  $$\begin{aligned}&{[}\text {pressure}_{N_\text {B31/B32},\text {ingestion}}] = [P_g > 0.5 p_\text {max}] \end{aligned}$$

  (53)
  $$\begin{aligned}&{[}\text {pressure}_{N_\text {B31/B32},\text {rejection}}] = [P_g > 0.25 p_\text {max}] \end{aligned}$$

  (54)
  $$\begin{aligned}&{[}\text {retracted}_{N_\text {B31/B32},\text {swallow},\text {off}}] = [x_{\text {g/h}}< z_{N_\text {B31/B32},\text {swallow},\text {off}}]\nonumber \\ \end{aligned}$$

  (55)
  $$\begin{aligned}&{[}\text {retracted}_{N_\text {B31/B32},\text {swallow},\text {on}}] = [x_{\text {g/h}}< z_{N_\text {B31/B32},\text {swallow},\text {on}}] \nonumber \\ \end{aligned}$$

  (56)
  $$\begin{aligned}&{[}\text {retracted}_{N_\text {B31/B32},\text {bite},\text {off}}] = [x_{\text {g/h}}< z_{N_\text {B31/B32},\text {bite},\text {off}}] \nonumber \\ \end{aligned}$$

  (57)
  $$\begin{aligned}&{[}\text {retracted}_{N_\text {B31/B32},\text {bite},\text {on}}] = [x_{\text {g/h}}< z_{N_\text {B31/B32},\text {bite},\text {on}}]\nonumber \\ \end{aligned}$$

  (58)
  $$\begin{aligned}&{[}\text {retracted}_{N_\text {B31/B32},\text {reject},\text {off}}] = [x_{\text {g/h}}< z_{N_\text {B31/B32},\text {reject},\text {off}}]\nonumber \\ \end{aligned}$$

  (59)
  $$\begin{aligned}&{[}\text {retracted}_{N_\text {B31/B32},\text {reject},\text {on}}] = [x_{\text {g/h}}< z_{N_\text {B31/B32},\text {reject},\text {on}}] \end{aligned}$$

  (60)

• 2. B6/B9/B3

  $$\begin{aligned} \begin{aligned}&N_\text {B6/B9/B3}(j+1) \\&\quad = N_\text {MCC}\,\,N_\text {B64}\,\,(\,!\,(N_\text {B4/B5}\ge 2)) \,\,\\&\bigg (\big ((N_\text {CBI-3}\,\,(\,!\,[\text {grasper}_\text {mech}])) \,\,[\text {pressure}_{N_\text {B6/B9/B3},\text {bite}}] \big ) + \\&\big ((N_\text {CBI-3}\,\,[\text {grasper}_\text {mech}]) \,\,[\text {pressure}_{N_\text {B6/B9/B3},\text {swallow}}] \big ) + \\&(\,!\,N_\text {CBI-3}) \,\,(\,!\,[\text {pressure}_{N_\text {B6/B9/B3},\text {reject}}]) \big ) \bigg ) \end{aligned} \end{aligned}$$

  (61)

  where

  $$\begin{aligned}&{[}\text {pressure}_{N_\text {B6/B9/B3},\text {bite}}] \nonumber \\&\quad = [P_g > z_{N_\text {B6/B9/B3},\text {bite},\text {pressure}})] \end{aligned}$$

  (62)
  $$\begin{aligned}&[\text {pressure}_{N_\text {B6/B9/B3},\text {swallow}}] \nonumber \\&\quad = [P_g > z_{N_\text {B6/B9/B3},\text {swallow},\text {pressure}})] \end{aligned}$$

  (63)
  $$\begin{aligned}&{[}\text {pressure}_{N_\text {B6/B9/B3},\text {reject}}] \nonumber \\&\quad = [P_g > z_{N_\text {B6/B9/B3},\text {reject},\text {pressure}}] \end{aligned}$$

  (64)

• 3. B8a/b

  $N_\text {B8}$ receives fast inhibitory and slow excitatory input from $N_\text {B40/B30}$ [33, 79]. In the Boolean framework here we implement this as an excitatory input immediately following cessation of $N_\text {B40/B30}$ activity for a user-specified duration ($\text {duration}_{N_\text {B40/B30},\text {excite}}$). Prior to calculating a new value for $N_\text {B8}$, we first check whether the synaptic connection from $N_\text {B40/B30}$ is excitatory with the following statements:

  if $(N_\text {B40/B30}(j) == 0 \, \text {AND} \, j < (t_{N_\text {B40/B30}\text {,off}} + \text {duration}_{N_\text {B40/B30},\text {excite}}))$, then set $N_\text {B40/B30},\text {excite} = 1$

  else set $N_\text {B40/B30},\text {excite} = 0$

  $$\begin{aligned} \begin{aligned} N_\text {B8}(j+1) =&N_\text {MCC}\,\,(\,!\,(N_\text {B4/B5}\ge 2)) \,\,\\&((N_\text {CBI-3}\,\,(N_\text {B20}\parallel (N_\text {B40/B30},\text {excite})) \,\,\\&\quad \quad (\,!\,N_\text {B31/B32})) + \\&\quad ((\,!\,N_\text {CBI-3}) \,\,N_\text {B20})) \end{aligned} \end{aligned}$$

  (65)

• 4. B7

  $$\begin{aligned} \begin{aligned}&N_\text {B7}(j+1) = N_\text {MCC}\\&\quad \,\,\big (\big ((\,!\,N_\text {CBI-3}\parallel [\text {grasper}_\text {mech}]) \\&\quad \,\,([\text {protracted}_{N_\text {B7},\text {reject}}] \parallel [\text {pressure}_{N_\text {B7}}])\big ) + \\&\big ((N_\text {CBI-3}\,\,\,!\,[\text {grasper}_\text {mech}]) \,\,([\text {protracted}_{N_\text {B7},\text {bite}}] \\&\quad \parallel [\text {pressure}_{N_\text {B7}}])\big )\big ) \end{aligned} \end{aligned}$$

  (66)

  where

  

  (67)
  

  (68)
  

  (69)

• 5. B38

  $$\begin{aligned} \begin{aligned} N_\text {B38}(j+1) = N_\text {MCC}\,\,[\text {grasper}_\text {mech}]\,\,\bigg ( N_\text {CBI-3}\,\,[\text {retracted}_{N_\text {B38}}] \bigg ) \end{aligned} \end{aligned}$$

  (70)

  where

  $$\begin{aligned} {[}\text {retracted}_{N_\text {B38}}] = [x_{\text {g/h}}< z_{\text {B38}}] \end{aligned}$$

  (71)

1.4 Muscle forces

Contact forces, such as the pressure resulting from grasper closure and force due to the anterior pinch, are implemented as second-order responses to neural activation using the semi-implicit integration scheme, Eq. (19), as shown in the following equations.

• 1. Grasper Pressure

  $$\begin{aligned}&P_\text {I4}(t+h) = \frac{\tau _\text {I4} P_\text {I4}(t) + h A_\text {I4}(t)}{\tau _\text {I4} + h} \end{aligned}$$

  (72)
  $$\begin{aligned}&A_\text {I4}(t+h) = \frac{\tau _\text {I4} A_\text {I4}(t) + h N_\text {B8}(t)}{\tau _\text {I4} + h} \end{aligned}$$

  (73)

• 2. Pinch Pressure

  $$\begin{aligned}&P_\text {I3,ant.}(t+h) = \frac{\tau _\text {I3,ant.} P_\text {I3,ant.} (t) + h A_\text {I3,ant.}(t) }{\tau _\text {I3,ant.} +h} \end{aligned}$$

  (74)
  $$\begin{aligned}&A_\text {I3,ant.}(t+h) \nonumber \\&\quad = \frac{\tau _\text {I3,ant.} A_\text {I3,ant.}(t) + h (N_\text {B38}(t) + N_\text {B6/B9/B3}(t))}{\tau _\text {I3,ant.} + h}\nonumber \\ \end{aligned}$$

  (75)

Muscle tensions for the remaining musculature were calculated using a second-order response to the neural activity as outlined in the following equations.

• 1. I3 Tension

  $$\begin{aligned}&T_\text {I3}(t+h) = \frac{\tau _\text {I3} T_\text {I3}(t) + h A_\text {I3}(t)}{\tau _\text {I3} + h} \end{aligned}$$

  (76)
  $$\begin{aligned}&A_\text {I3}(t+h) = \frac{\tau _\text {I3} A_\text {I3}(t) + h N_\text {B6/B9/B3}(t)}{\tau _\text {I3} + h} \end{aligned}$$

  (77)

• 2. I2 Tension

  Time constants for I2 were tuned independently for ingestion and egestion to account for the experimental observations that egestions have a longer period than ingestions. Such variation in responsiveness of the animal may exist due to differences in neuromodulation between the behaviors. Therefore, the time constant for I2 is calculated as:

  $$\begin{aligned}&\tau _\text {I2} =N_\text {CBI-3}\;\tau _\text {I2,ingestion} \nonumber \\&\qquad + (\,!\,N_\text {CBI-3}) \;\tau _\text {I2,egestion} \end{aligned}$$

  (78)
  $$\begin{aligned}&T_\text {I2}(t+h) = \frac{\tau _\text {I2} T_\text {I2} + h A_\text {I2}}{\tau _\text {I2} + h} \end{aligned}$$

  (79)
  $$\begin{aligned}&A_\text {I2}(t+h) = \frac{\tau _\text {I2} A_\text {I2}(t) + h N_\text {B31/B32}(t)}{\tau _\text {I2} + h} \end{aligned}$$

  (80)

• 3. Hinge Tension

  $$\begin{aligned} T_\text {hinge}(t+h)= & {} \frac{\tau _\text {hinge} T_\text {hinge}(t) + h A_\text {hinge}(t)}{\tau _\text {hinge} + h} \end{aligned}$$

  (81)
  $$\begin{aligned} A_\text {hinge}(t+h)= & {} \frac{\tau _\text {hinge} A_\text {hinge}(t) + h N_\text {B7}(t)}{\tau _\text {hinge} + h} \end{aligned}$$

  (82)

1.5 Biomechanical model

The motions of the head and grasper are calculated based on the quasi-static equations of motion:

$$\begin{aligned} \frac{\mathrm{d}}{\mathrm{d}t} \begin{bmatrix} x_\text {g}\\ x_\text {h}\end{bmatrix}= {\begin{bmatrix} \frac{F_g}{c_g}\\ \frac{F_h}{c_h} \end{bmatrix}} \end{aligned}$$

(83)

where $x_\text {h}$ is the position of the head relative to the ground frame, $x_\text {g}$ is the position of the grasper relative to the ground frame, and $c_h$ and $c_g$ are the damping coefficients for the motion of the head and grasper, respectively. The forces on the grasper and head can be calculated as outlined in the following sections.

1.5.1 Forces on the grasper

The positive direction for $x_\text {g}$ corresponds to protraction (Fig. 3). The sum of forces on the grasper is

$$\begin{aligned} F_g = {F}_\text {I2} + F_{sp,g} - {F}_\text {I3} - {F}_\text {hinge} + {F}_{f,g} \end{aligned}$$

(84)

where the component forces are defined and calculated as follows:

${F}_\text {I2}$: The force due to the I2 muscle. This value is dependent on the tension of the muscle as well as the mechanical advantage. It is scaled by a tunable maximum parameter, $F_\text {I2,max}$, and is calculated as follows:

$$\begin{aligned} {F}_\text {I2} = F_\text {I2,max} T_\text {I2}(t)(1-x_{\text {g/h}}) \end{aligned}$$

(85)

where $x_{\text {g/h}}= x_\text {g}- x_\text {h}$ is the position of the grasper relative to the head.

$F_{sp,g}$: The force in the spring connecting the grasper to the head. This spring represents the surrounding musculature of the esophagus, buccal mass, and extrinsic muscles which are not explicitly modeled. This is calculated as:

$$\begin{aligned} F_{sp,g} = K_g (x_{g/h}^0-x_{\text {g/h}}) \end{aligned}$$

(86)

where $K_g$ is the spring constant and $x_{g/h}^0$ is the rest length of the spring.

${F}_\text {I3}$: The force due to the I3 muscle which pushes the grasper backwards during retraction. This force is due to tension in I3 closing the muscular toroids. This value is dependent on the tension of the muscle as well as the mechanical advantage. It is scaled by a tunable maximum parameter, $F_\text {I3,max}$, and is calculated as follows:

$$\begin{aligned} {F}_\text {I3} = F_\text {I3,max} T_\text {I3}(t)(x_{\text {g/h}}-0) \end{aligned}$$

(87)

${F}_\text {hinge}$: The force due to the hinge. This value is dependent on the tension of the muscle as well as the mechanical advantage. It is scaled by a tunable maximum parameter, $F_\text {hinge,max}$, and is calculated as follows:

$$\begin{aligned} {F}_\text {hinge} = [\text {hinge stretched}]F_\text {hinge,max} T_\text {hinge}(t)(x_{\text {g/h}}-0.5) \end{aligned}$$

(88)

where $[\text {hinge stretched}] = [x_{\text {g/h}}>0.5]$ determines whether the hinge is sufficiently stretched to produce any force [148].

$F_{f,g}$: Friction resulting from the grasper closing on an object. To determine $F_{f,g}$ it is necessary to check if the grasper is slipping against the object by checking the inequality:

$$\begin{aligned} |{F}_\text {I2} + F_{sp,g} - {F}_\text {I3} - {F}_\text {hinge}| \le |\mu _{s,g} {F}_\text {I4}| \end{aligned}$$

(89)

where $\mu _{s,g}$ is the coefficient for static friction between the grasper and the object. ${F}_\text {I4}$ is the normal force due to the grasper muscle I4 closing on the object. This is calculated directly as the grasper pressure defined in the previous appendix applied to a unit area scaled by a parameter.

$$\begin{aligned} {F}_\text {I4} = F_{I4,\text {max}} P_\text {I4}(t) \end{aligned}$$

(90)

If the condition in Eq. (89) is true, then the contact is in a state of static friction and $F_{f,g}$ is calculated as:

$$\begin{aligned} |F_{f,g}| = {F}_\text {I2} + F_{sp,g} - {F}_\text {I3} - {F}_\text {hinge} \end{aligned}$$

(91)

If the condition in Eq. (89) is not true, the contact is sliding and is in a state of kinetic friction, and $F_{f,g}$ is calculated as:

$$\begin{aligned} |F_{f,g}| = \mu _{k,g} {F}_\text {I4} \end{aligned}$$

(92)

where $\mu _{k,g}$ is the coefficient for kinetic friction between the grasper and the seaweed.

The sign of the friction force is dependent on which direction the grasper would be moving without the friction present, and ${F}_{f,g}$ can be calculated as:

$$\begin{aligned} {F}_{f,g} = -\text {sgn}({F}_\text {I2} + F_{sp,g} - {F}_\text {I3} - {F}_\text {hinge}) |F_{f,g}| \end{aligned}$$

(93)

1.5.2 Forces on head

The forces on the head are calculated as:

$$\begin{aligned} F_h = F_{sp,h} - F_{sp,g} - F_\mathrm{I2} + F_\mathrm{I3} + F_\mathrm{hinge} + F_{f,h} \end{aligned}$$

(94)

The muscles and grasper spring exert forces on the head equal and opposite to those on the grasper. As the muscles contract and apply forces to move the grasper forward, this also stretches the spring between the grasper and head proportionally to the muscle force. For the quasi-static model, acceleration is assumed to be negligible and therefore the forces on the grasper must equal zero.

$$\begin{aligned} 0 = F_{sp,g} + F_\text {I2} - F_\text {I3} - F_\text {hinge} + F_{f,g} \end{aligned}$$

(95)

Solving for the spring forces, $F_{sp,g}$, and substituting into Eq. (94) yields:

(96)

which simplifies to:

$$\begin{aligned} F_h = F_{sp,h} + F_{f,g} + F_{f,h} \end{aligned}$$

(97)

where $F_{sp,h}$ is the spring force between the head and neck of the animal, $F_{f,g}$ is the previously calculated friction force between the grasper and the object, and $F_{f,h}$ is the friction force resulting from the jaws pinching on the object. These components are calculated as follows.

$$\begin{aligned} F_{sp,h} = K_h (x_{h}^0-x_{h}) \end{aligned}$$

(98)

where $K_h$ is the spring constant and $x_{h}^0$ is the rest length of the spring.

To determine the value of $F_{f,h}$ it is necessary to check if the jaws are slipping relative to the seaweed by checking the following inequality:

$$\begin{aligned} |F_{sp,h} + F_{f,g}| \le |\mu _{s,h} {F}_\text {I3,ant.}| \end{aligned}$$

(99)

where $\mu _{s,h}$ is the coefficient for static friction between the jaws and the seaweed. ${F}_\text {I3,ant.}$ is the normal force due to the anterior portion of the I3 jaw muscle closing on the seaweed. This is calculated directly as the pinch pressure defined in the previous appendix applied to a unit area, scaled by a parameter, and multiplied by a mechanical advantage term:

$$\begin{aligned} {F}_\text {I3,ant.} = F_{\text {I3,ant.,max}} P_\text {I3,ant.}(t)(1-x_{\text {g/h}}). \end{aligned}$$

(100)

If the condition in Eq. (99) is true, the jaws are in static friction and $F_{f,h}$ is calculated as:

$$\begin{aligned} |F_{f,h}| = F_{sp,h} + F_{f,g} \end{aligned}$$

(101)

If the condition in Eq. (99) is not true, the jaws are slipping and $F_{f,h}$ is calculated as:

$$\begin{aligned} |F_{f,h}| = \mu _{k,h} {F}_\text {I3,ant.} \end{aligned}$$

(102)

where $\mu _{k,h}$ is the coefficient for kinetic friction between the jaws and the seaweed.

The sign of the friction force is dependent on which direction the head would be moving without the friction present and ${F}_{f,h}$ can be calculated as:

$$\begin{aligned} {F}_{f,h} = -\text {sgn}(F_{sp,h} + F_{f,g}) |F_{f,h}| \end{aligned}$$

(103)

1.5.3 Force on objects

The force on the object if unbroken is equal to the sum of the friction forces where we use the conventions that positive force indicates tension on the force transducer:

$$\begin{aligned} F_o = F_{f,g} + F_{f,h} \end{aligned}$$

(104)

If $F_o \le z_s$, where $z_s$ is the user defined seaweed strength, the seaweed is not broken and the motion of the bodies is calculated based on the forces calculated in the previous sections. If $F_o > z_s$, the seaweed is broken and can no longer transmit forces to the head or grasper. Therefore, the forces on the head and grasper are recalculated as:

$$\begin{aligned} F_h= & {} F_{sp,h} \end{aligned}$$

(105)

$$\begin{aligned} F_g= & {} {F}_\text {I2} + F_{sp,g} - {F}_\text {I3} - {F}_\text {hinge} \end{aligned}$$

(106)

A Boolean tracking variable [unbroken] is used to track whether the seaweed is intact (1) or broken (0). Once the seaweed breaks, it is not restored until the grasper has completed a new protraction and grasp motion. For this model, we have implemented this by resetting $[\text {unbroken}] = 1$ if at the current timestep $[\text {unbroken}] ==0$ AND $x_{\text {g/h}}< 0.3$ AND $x_{\text {g/h}}(j+1)>x_{\text {g/h}}(j)$. These thresholds were tuned manually for this implementation.

1.5.4 Updating grasper and head positions

All of the forces in this biomechanical model are linearly dependent on the position of the head, $x_\text {h}$, and grasper, $x_\text {g}$. Therefore, they can each be rewritten in the form:

$$\begin{aligned} F = \mathbf{A} _F \begin{bmatrix} x_\text {g}\\ x_\text {h}\end{bmatrix}+ \mathbf{b} _F \end{aligned}$$

(107)

As a consequence, the equations of motion can be rewritten in the form

$$\begin{aligned} \frac{\mathrm{d}}{\mathrm{d}t} \begin{bmatrix} x_\text {g}\\ x_\text {h}\end{bmatrix}= {\begin{bmatrix} \frac{A_{11}}{c_h} &{} \frac{A_{12}}{c_h}\\ \frac{A_{21}}{c_g} &{} \frac{A_{22}}{c_g} \end{bmatrix}} \begin{bmatrix} x_\text {g}\\ x_\text {h}\end{bmatrix}+ \begin{bmatrix} b_{1} \\ b_{2} \end{bmatrix} \end{aligned}$$

(108)

This can then be integrated with the semi-implicit integration scheme in Appendix A.1 as:

$$\begin{aligned}&\mathbf {x}(t+h) = \frac{1}{1-h\text {Tr}A(t)}\nonumber \\&\quad \left[ \left( I+h\left( \begin{array}{rr} -A_{22} &{} A_{12}\\ A_{21} &{} -A_{11} \end{array} \right) \right) \mathbf {x}(t)+h\mathbf {b}(t)\right] . \end{aligned}$$

(109)