In vivo electrophysiology recordings and computational modeling can predict octopus arm movement

Abstract

The octopus has many features that make it advantageous for revealing principles of motor circuits and control and predicting behavior. Here, an array of carbon electrodes providing single-unit electrophysiology recordings were implanted into the octopus anterior nerve cord. The number of spikes and arm movements in response to stimulation at different locations along the arm were recorded. We observed that the number of spikes occurring within the first 100 ms after stimulation were predictive of the resultant movement response. Machine learning models showed that temporal electrophysiological features could be used to predict whether an arm movement occurred with 88.64% confidence, and if it was a lateral arm movement or a grasping motion with 75.45% confidence. Both supervised and unsupervised methods were applied to gain streaming measurements of octopus arm movements and how their motor circuitry produces rich movement types in real time. For kinematic analysis, deep learning models and unsupervised dimensionality reduction identified a consistent set of features that could be used to distinguish different types of arm movements. The neural circuits and the computational models identified here generated predictions for how to evoke a particular, complex movement in an orchestrated sequence for an individual motor circuit. This study demonstrates how real-time motor behaviors can be predicted and distinguished, contributing to the development of brain-machine interfaces. The ability to accurately model and predict complex movement patterns has broad implications for advancing technologies in robotics, neuroprosthetics, and artificial intelligence, paving the way for more sophisticated and adaptable systems.

Introduction

Current Brain Machine Interface (BMI) systems are based on decoding algorithms that use the neural signals to control the external device (Carmena et al. 2003; Downey et al. 2016). However, these devices do not provide enough independent degrees of freedom of the arm, and usually control even simple motions of the lower limbs. Understanding the motor circuits that govern the trajectory and dynamics of arm movement is crucial for advancing neuroprosthetic devices and robotic limbs capable of precise actions like reaching and grasping. Recent studies using different computational techniques including machine learning (ML) algorithms and Artificial Intelligence (AI) have shown the ability to predict several aspects of arm reaching from electrophysiology data (Yoshimura et al. 2017; Pirondini et al. 2017; Kumarasinghe et al. 2021). To improve future BMI devices, it will be crucial to further reveal the neural mechanisms behind how diverse movements are represented in the measured electrophysiological signals and how these representations relate to distinct kinematic features of the behavioral response (position, velocity, muscle activity, direction, and more).

It remains challenging to identify specific motor control computational circuits because the nervous system has a complex and dense neural architecture. In addition, some cognitive neuroscience theories suggest that motor decisions incorporates computational circuits outside the central nervous system (Bluhm 2017; Informatics and the Search for Mental Structure 2016; RA, P. 2010). In order to build predictive models it will be required to record neural correlates of movement from the brain as well as from local peripheral circuits. Studying motor control in species with a decentralized nervous system can reveal fundamental and alternative computational paradigm of motor control.

The octopus has many features that makes it advantageous for pursuing a holistic understanding of movement, cognition and behavior: It has a highly developed nervous system containing 500 million neurons and a large brain (Wells 1978; Young 1971; VanBuren et al. 2021); Each of the eight arms contains an axial nerve cord (ANC) which resembles and acts like the vertebrate’s spinal cord; Hundreds of suckers and sucker ganglia act as the peripheral nervous system and can demonstrate a large repertoire of behaviors (Jarmoluk and Pelled 2024; Mather and Alupay 2016); it has distributed control of its nervous system. One of the advantages of using this animal model, is that an electric or tactile stimulation to the octopus’s denervated arm, can still trigger movement that is similar in kinematics to movement of an intact arm as was previously demonstrated (Richter et al. 2015; Sumbre et al. 2001, 2006; Hanassy et al. 2015). This suggests that the octopus has a simplified neural program embedded within the arm itself and is adaptable to various degrees of input from visual, sensory and motor brain areas (Richter et al. 2015; Zullo et al. 2009, 2019). Classifications of movement in mammals (e.g., rats, humans) have revealed fundamental principles of sensorimotor control and coordination. The octopus, with its decentralized nervous system and relatively simple neural circuits, provides a unique model for studying these principles. Unlike mammals, whose complex and overlapping neural networks can make it challenging to isolate specific pathways, the octopus’s less centralized nervous system allows for a clearer understanding of core adaptive strategies and motor pathways. This simplicity presents an opportunity to apply machine learning and advanced statistical methods to develop computational models capable of predicting motor behavior in a way that may be more challenging in mammals.

Recently, there has been growing interest in recording electrophysiology signals from octopus nervous system: single spike and intracellular recordings from slices obtained from octopus had revealed principles of learning and memory (Nesher et al. 2019), and multi-unit and local-filed potentials were recorded from octopus arms (Zullo et al. 2019) and the nerve ring responsible for arm coordination (Chang and Hale 2023). Recent developments have also shown the capabilities of recording brain signals from awake octopuses (Gutnick, et al. 2023). These studies revealed fundamental concepts regarding the octopus nervous system, motor control and coordinated movement. However, thus far, the electrophysiology recordings have only consisted of continuous local field and electroencephalogram data. For this study, as described below, we instead used very small diameter carbon fiber electrodes for spike recording. Carbon fiber electrodes are strong enough to penetrate tough neural tissue when sharpened (Welle et al. 2021; Letner, et al. 2023). They also do very little damage to the neurons of interest (Patel et al. 2016). They can be sharpened in a way that preserves a small electrode surface area, enabling high amplitude spikes (Richie, et al. 2024).

While the octopus demonstrates useful complex and flexible movements, this kinematics must first be measured and quantified for correlation with electrophysiological signals. Analyzing movements in an automated manner can quickly provide crucial insights into animal behavior that would otherwise be too time-intensive or costly to manually characterize. In many situations, however, a pre-defined quantity such as spatial location may not be the best metric for characterizing behavior. More complex and ethologically relevant behaviors, such as exploration, reaching, or grasping, are better defined by their motion with respect to the animal’s body or to the sequence in which they are performed (McCullough and Goodhill 2021). Thus, methods are needed for automated position extraction, pose estimation, and behavioral feature identification.

Both supervised and unsupervised machine learning methods can be used across a wide variety of animal models to classify and cluster behavioral features into these more relevant phenotypes (Wiltschko and Wiltschko 2005; Marques et al. 2018; Hsu and Yttri 2021). Supervised machine learning effectively quantifies behavioral data, while unsupervised clustering objectively uncovers inherent structures within datasets, aiding in identifying continuous movement spaces or distinct movement type clusters. While characterizing behavior itself is important for understanding animals and their nervous systems, it lacks the ability to test hypothesized brain-behavior links. That analysis requires models that are able to correlate neural activity and behavioral responses to perturbations. By learning behavioral features associated with various experimental paradigms, we could then correlate what aspects of the environment or stimuli are significantly driving these behaviors. Identifying a set of movements and their orchestrated sequences empowers the construction of simplified yet accurate representations for a particular task, shedding light on underlying mechanisms of e.g., the motor circuits involved in reaching.Finally, with the development of tools that allow real-time analysis with minimal latency (Mathis et al. 2018; Draelos et al. 2021; Kane et al. 2020), we can also consider closed-loop experimental paradigms that adapt stimulation parameters based on instantaneous behavioral responses. With immediate analysis of how different behavioral features vary during neural stimulation, we could construct models that learn how best to evoke a particular, complex movement in an orchestrated sequence for a particular motor circuit.

This type of data-driven approach could unveil individual behavioral motifs, control circuits, and ultimately contribute to advancements in more flexible and adaptable prosthetics; notably in goal-oriented grasping movements for individuals with limb loss or spinal cord injuries. Here, using state-of-the-art carbon fiber arrays that provided single-unit and multi-unit electrophysiology recording capabilities of ANC neurons (Richie, et al. 2024), we obtained simultaneous video recordings, neural patterns and arm kinematics. To trigger movement, descending stimulation was delivered directly on the ANC, and ascending tactile stimulation was delivered to the base of the arm, close to the electrodes, and to the more distal portion of the arm. Machine learning models were built to then predict resultant octopus arm movement and learn what specific behaviors could be decoded to the original stimulation.

Methods

Experiments & data acquisition

All procedures were approved by the Institutional Animal Care and Use Committee at Michigan State University. Adult Octopus bimaculoides (n = 7) were anesthetized according to procedures described by Abbo et. al., (2021). Ethanol was chosen as the anesthetic agent over magnesium chloride due to its shorter induction and recovery times. The octopuses were anesthetized in a 6L plastic container that contained 4L of water from their home tanks and was aerated using an air stone. 50 mL of 200 proof pure ethanol was added to the water just before the octopus was placed in the container for a starting concentration of 1.25% ethanol. Breath rate was used as the primary means of determining when the anesthesia process was complete, along with additional indicators as described by Fiorito et al., (2015). Additional ethanol was added in 5-min increments as needed for a total concentration of no more than 3% ethanol until the octopus’s breathing rate fell between 10–20 breaths/minute and the animal was nonreactive to noxious stimuli in both its central (mantle pinch) and peripheral (distal arm pinch) nervous systems. From each animal we recorded from two arms (L2 and R3), and the recordings of each of the arms was spaced at least three weeks apart as only one arm was removed at a time during each procedure. Once the arm was removed, and the proximal end of the arm was restrained in a tray that was continuously perfused with filtered saltwater. The suckers had ventral orientation and the stimulations occurred on the dorsal side of the arm. The muscles at the base of the arm were dissected, revealing the ANC. In house high-density carbon fiber of 16 electrodes array was inserted transversely into the exposed ANC (Welle et al. 2021; Richie et al. 2024), targeting the cerebrobrachial tracts to induce efferent signaling (Fig. 1). To determine if there was a difference in stimulation response, mechanical (tactile) stimulation was performed using plastic forceps and an applied force of 2-3lbf in the same location as the electrical stimulations. Electrical stimulation was performed using a single bipolar electrode driven by an IsostimA320 that delivered 120v 5 mA, 5 pulses with a duration of 5 ms. The arm was stimulated in three different locations: directly on the exposed ANC to reflect efferent stimulation, and on regions proximal and distal to the electrode placement, reflecting afferent stimulation. A 2-sample Kolmogorov–Smirnov test was used to compare the distributions of kinematic features between the tactile and electrical stimulation conditions. The results indicated that we could not reject the hypothesis that the distributions were drawn from the same underlying population. Consequently, we concluded that the features we examined did not provide a statistically significant basis for distinguishing between the two stimulation types. Therefore, all stimulation trials in each of the locations were grouped.

Fig. 1

The experiment setup. A The 16-electrode array was inserted into the exposed axial nerve cord (ANC). Tactile and electrical stimulation were delivered into three different locations: directly onto the ANC, proximal, and distal to the electrodes’ placements. The carbon electrode array recorded spikes. B Representative spikes are shown on the right

Intan Recording System (Intan Technologies, Los Angeles, CA) and Spike2 (Cambridge Electronic Design, Cambridge, UK) were used to record signals. Spike2 was used for initial signal processing. The signals were filtered with a bandpass second order Butterworth 0.1 to 3 kHz and the threshold for action potential (AP) detection was set at 4 standard deviations (Welle et al. 2020; Lu et al. 2008; Pelled et al. 2009). APs had an average duration of 1.44 ms. A total of 95 experiments resulted in 1520 traces of recordings. For modeling analysis, recordings from three electrodes in the array that showed the highest activity were selected (275 traces).

For movement recording, a webcam was positioned over the recording chamber and videos were recorded simultaneously with electrophysiology recordings via the Spike2 interface. Movement was first manually classified into 3 distinct categories: no movement (NM, “0”), lateral movement (LM, “1”), and a curl (CM, “2”).

Electrophysiology data of the entire 1520-spike recordings traces was processed using Plexon (Plexon Inc, Dallas, TX) offline sorting software with 250 Hz Butterworth filters applied to the raw data. Spikes that had passed the threshold were identified as “spikes” in the analysis. The processed data was analyzed using Python. The Python script extracted the number of spikes for each 50 ms time bin of the 16 channels and summed them. Prism software (Graphpad software Inc, San Diego, CA) was used for statistics. Outlier data was filtered using the 1% ROUT method. Welch and Brown-Forsythe ANOVA test was used to assess the statistical significance of the relationships between spikes and movement, spikes and stimulation location, and location of stimulation and type of movement. Statistically significant results were considered to be p < 0.05.

Modeling to predict movement from neural activity

We utilized the One Hot Encoding (OHE) technique (Yang et al. 2018; Al-Shehari and Alsowail 2021) to convert the categorical feature into binary features. The OHE technique converts a categorical feature into multiple binary features, where the number of binary features is equivalent to the number of distinct categories in the original categorical feature. This technique assigns a value of 1 to the binary feature corresponding to the specific category for each instance/sample and a value of 0 to all other binary features.

We utilized Cramer's V (Cramér 1946) to understand the significance of each OHE feature with the categorical target. Cramer's V is a technique used to measure the degree of association between two categorical features. This technique is based on the chi-square statistic test, and Cramer's V value ranges from 0 to 1, where 0 indicates no association between the variables, and 1 indicates a perfect association between the variables.

An additional machine learning based method, the Feature importance analysis (Xu et al. 2023; Moslehi et al. 2022) can reveal the degree of importance of all features, including categorical and discrete, binned, features to predict the target. The analysis was conducted on both Binary-class (movement/no movement) and Multi-class (no movement/movement/movement with a curl) datasets to identify which features were most influential in predicting the movement outcomes. This analysis was essential to understand the underlying factors contributing to the model's predictions, optimize the model's performance, and provide insights into the key drivers of movement patterns.

Overall, an array of 16 different machine learning models was trained on our rich datasets, demonstrating a comprehensive application of diverse machine learning techniques across several categories. Tree-based models like the Decision Tree, ensemble techniques such as robust methods like Random Forest and Extra Trees Classifier, as well as powerful boosting approaches such as Extreme Gradient Boosting (XGBoost), Light Gradient Boosting Machine (LightGBM), Gradient Boosting Classifier, Adaptive Boosting (Ada Boost), and the state-of-the-art CatBoost were employed. Advanced classifiers like the Ridge Classifier, which is based on linear regression techniques but includes regularization, and SVM with a Linear Kernel, which excels in high-dimensional spaces, were also employed. Additionally, simpler yet essential models like the Dummy Classifier were used to establish baseline performances. Our model set further incorporated classical statistical methods, including Logistic Regression and Naive Bayes, alongside discriminant analysis techniques like Linear and Quadratic Discriminant Analysis. The array also included instance-based learning methods like K-Nearest Neighbors (KNN). This varied and methodologically rich collection of models ensured rigorous and nuanced analysis, providing robust and detailed insights.

Modeling to identify stimulation from resultant behavior

To track the motion of the octopus arm, we first employed DeepLabCut (DLC) for markerless keypoint tracking and pose estimation (Mathis et al. 2018; Kane et al. 2020). This widely-used software package utilizes deep neural networks and transfer learning to achieve accurate 2D and 3D markerless pose estimation for defining and tracking specific points of interest. Out of the total 234 videos, 16 videos with different camera angles, applied stimuli, and observed motion types were selected to train our octopus-specific model. 16 images from each of these videos were then selected by DLC as representative and diverse samples of the octopus arm’s movement, as determined by k-means clustering. Images that were blurred and where the octopus was heavily obscured were then manually dropped from the training set (typically 0–4 images per video). Finally, the images were hand-annotated to label 17 (approximately) equidistant key points along the arm using a GUI provided by the DLC package (Fig. 1A).

Next, we took the ResNet-50 model supplied by DLC, pre-trained on the large and well-established ImageNet dataset, and further trained it using our annotated frames of the octopus arm. This training was done on a lab workstation with a single GPU (NVIDIA GeForce RTX 4070 Ti) and took 3.5 h to run 150,000 iterations. Once the training was complete, the final model was employed for real-time key point tracking of all videos in the dataset. The model output reported x and y location predictions for each of the 17 key points in each frame, accompanied by an associated prediction likelihood value.

To comprehensively quantify the entire arm’s motion, a range of significant kinematic features were computed from the x and y predictions. Specifically, we defined θ as the angle formed by the proximal and distal segments (the angle between the stationary base and the tip of the arm), and its instantaneous angular speed as a difference of θ across consecutive frames scaled by frame-per-second to convert to SI spikes. We also considered the absolute speed of each keypoint, later focusing on just the distal point as significant. To provide an initial quantification across time, the mean and maximum values of the above features were calculated over three distinct non-overlapping time intervals after stimulation: 0–1 s, 1–2 s, and 2 or more seconds (Fig. 2b). These defined features can be well understood, linked as they are to specific locations along the arm. For example, the distribution of the maximum angular velocity of the arm in the 0–1 s time period post-stimulation is clearly different for the videos hand-labeled as having ‘No movement’, ‘Movement’, and ‘Movement with arm curl’ (Fig. 2c), as might be intuitively expected.

Fig. 2

Keypoint tracking and feature extraction. a 17 roughly equidistant points were hand-labeled along the length of the arm. Using these keypoint positions, various metrics were computed for subsequent analysis including the overall angle made between the (stationary) base and the tip (θ) and the angular and keypoint velocities. The distal-most keypoint (x) and its velocity (v) were found to be significant in distinguishing motion types. b To quantify motion across time, three intervals post-stimulation were considered: the first (t0) and second (t1) seconds, where most motion occurred, and 2 or more seconds (t2) until any observed motion ceased. c The histogram of one example metric, the maximum angular velocity in t0, is plotted (using a kernel density estimator) for each of the human-labeled movement categories. As expected, the ‘No movement’ videos have very low or zero angular velocity, whereas ‘Movement with arm curl’ videos tend to have higher maximum angular velocity

We then employed a second method of analysis that did not rely on key points, as typical key point tools rely on obvious features such as joints or consistent markings that the octopus arm lacks. We used an unsupervised streaming dimension reduction algorithm known as Procrustean SVD (proSVD) (Draelos et al. 2021) to identify features within the videos and how they varied across time, without any pre-training or knowledge of what the videos contained. Unlike conventional SVD methods, proSVD stands out by ensuring the selection of a stable feature set across time, offering dependable results even in the initial phases of data collection. We reduced the videos to 4 bases, or features, and quantified the discovered motion with the L2-norm of each basis vector. Additionally, to optimize processing efficiency, specialized code was developed to crop the videos precisely around the identified DLC key points with a 20 pixel margin before handing the cropped videos to proSVD. This tailored step proved instrumental in eliminating superfluous background elements (anything not an octopus arm), which significantly sped up subsequent processing stages.

Results

Single and multiunit analysis

Carbon fiber electrodes successfully recorded single and multiunit activity from the ANC, as shown in Fig. 1. The total number of spikes in the first 50 ms and 100 ms after stimulation were calculated for each of the 16 channels, which resulted in 1520 traces. The movement was based on video analysis and was classified into 3 distinct responses: no movement (“NM”), lateral movement (“LM”), and a movement that consisted an arm curl (“CM”). First, we measured the number of spikes occurring between 0-50 ms and 50-100 ms during each movement. All the stimulation locations were grouped. Figure 3 shows the number of spikes for each movement. To test if the number of spikes occurring immediately after stimulation is different for each movement response, a Welch and Brown-Forsythe ANOVA analysis was performed. Results showed that there was a significant difference between the groups means for each movement response type (F(2.00, 42.09) = 4.10, p = 0.023)). There was also a significant difference between the number of spikes occurring 50-100 ms after stimulation (F(2.00, 46.37) = 7.36, p = 0.0017). In both time frames, the lateral movement showed the greatest number of spikes activity. This may suggest that an arm curl is a reflexive response which requires less neural activation.

Fig. 3 

The number of spikes and movement. A The average number of spike responses within 50 ms after stimulation evoked different movement response. B The average number of spike responses within 50-100 ms after stimulation and the evoked movement response. In both time periods, the greatest number of spikes was found to be associated with lateral movement (LM). (No movement (MN), movement with an arm curl (CM); Welch and Brown-Forsythe ANOVA, * < 0.05; ** < 0.005; n = 7 octopuses, data obtained form 14 arms)

We then tested if there was a difference between the number of spikes occurring as a response to the location of the stimulation. The results show in Fig. 4 demonstrate that in the first 50 ms and 100 ms after stimulation there is a significant difference between the groups means (F(2.00, 54.23) = 6.062, p = 0.0042) and (F(2.00, 68.75) = 4.72, p = 0.012), respectively. In both time frames, the Cord stimulation showed the least variance in the number of spikes, compared to Distal and Proximal arm stimulation, suggesting that afferent stimulation results in a consistent response. Results also show that the number of spikes evoked by a Distal stimulation significantly increases over time which may suggest a mechanism to amplify distant signals (paired T-test analysis, p = 0.0067). The number of spikes between the first 50 ms and 100 ms did not change in response to Proximal or Cord stimulation.

Fig. 4

The stimulation type and number of spikes. A The average number of spike responses within 50 ms after stimulation and the location of the stimulation. B The average number of spike responses within 50-100 ms after stimulation and the location of the stimulation. Between the two time periods, the greatest increase in spikes was found to be associated with distal stimulations. (Welch and Brown-Forsythe ANOVA, * < 0.05; ** < 0.005; n = 7 octopuses, data obtained form 14 arms)

We then sought to determine the probability of the type of stimulation to evoke a specific movement response. Figure 5 shows the probability of movement response given the type of stimulation. Distal stimulation showed a clear preference to induce movement; In 94% of trials, it evoked a lateral movement (41%) or an arm curl (53%). On the other hand, Proximal and Cord stimulations did not induce a consistent response: Proximal stimulation induced lateral movement (25%), arm curl (29%), and in 46% of trial no movement was evoked; Cord stimulation induced lateral movement (31%), arm curl (31%), and in 38% of trial no movement was evoked.

Fig. 5

The stimulation type and movement response. A The probability of evoking each movement to occur after a distal stimulation. B The probability of evoking each movement to occur after a proximal stimulation. C The probability of evoking each movement to occur after a Cord stimulation. The variability of the movement response within each arm is shown for each stimulation location. This reveals that distal stimulation was very likely to induce a movement, whereas cord and proximal had a more even distribution of responses. (No movement (MN), movement with an arm curl (CM), lateral movement (LM)

The results showed that the number of spikes in the first 100 ms post-stimulation can predict the movement response. The predictive probability of longer period of spikes to inform on movement response was examined. The total number of spikes in the first 500 ms after stimulation was calculated for each of the 16 channels. The average number of spikes in the first 500 ms for movement MN was 591±33.2, LM resulted in 691±40 spikes, and CM resulted in 706±54.7 spikes (average±SEM). A Welch and Brown-Forsythe ANOVA analysis showed that there wasn’t a significant difference between the number of spikes to the movement response (F(2.00, 18.56) = 0.39, p = 0.68). In addition, there wasn’t a significant difference between the stimulation type and the number of spikes (F(2.00, 64.22) = 0.45, p = 0.64).

Computational modelling of electrophysiology responses

Feature extraction

Temporal and stimulation features were extracted from the electrophysiology signals, using sample windows of 3 s long, binned into 100 ms, and summing the number of spikes in each window. The ML analysis looked for correlations as a function of time and included more complex and longer data sets. 3 s was chosen to determine which time bins were the most important for further analysis. This process resulted in a dataset with 30 bins that we treated as discrete features. To encode the stimulation information, categorical features were added using OHE technique to convert the categorical data into a numerical format. Two different datasets of 275 traces were created: a Multi-Class dataset where each sample was labelled 0 (no movement; 74 samples), 1 (movement; 96 samples), and 2 (movement with arm curl; 105 samples); and a Binary dataset where samples labelled 1 and 2 were combined into a single movement class (1, consisting of 201 samples), and the 0 class (74 samples). The distribution of the samples in the multi-class dataset was found to be balanced in the number of trials in each class. However, the distribution of the samples in the binary dataset showed a slight imbalance as it consisted of more samples in movement 1. A flow chart describing the ML pipeline is illustrated in Fig. 6.

Fig. 6

An illustration of the machine learning pipeline

Next, we created a dataset by extracting temporal and stimulation features from the electrophysiology signals which contained 34 features: 33 features were predictors, and one feature was the categorical movement response (Target). Among 33 predictors, 30 predictors were discrete features, and the remaining three predictors were OHE features derived from the stimulation location: ANC stimulation, and stimulations located in the distal part and proximal part of the arm.

A Cramer's V analysis was computed to understand the impact of binary features as shown in Table 1. The stimulation type feature was encoded using OHE technique. Then, Cramer's V was computed between each binary feature and the target feature, reflecting the degree of association. Significantly, for both binary and multi-class targets, the tactile, distal part binary feature is highly associated with the target, indicating that the tactile, distal part binary feature is more useful for predicting the target features. The relatively low values for the ANC binary feature suggest it might be less useful for predicting the target features.

Table 1 Association strengths between binary features and target outcomes using Cramer’s V

Mutual information is a statistical method measuring the amount of dependence between two random variables (Witten et al. 2016), and Fig. 7 depicts the scores for both binary-class and multi-class.

Fig. 7

Mutual information analysis for understanding feature importance in binary and multi-class classification. A Mutual Information Scores between 30 binned discrete features, and Binary-Class scenario and (B) Mutual Information Scores between 30 binned discrete features, and Multi-Class scenario. Color bar indicates the mutual information score for a feature. This analysis demonstrates that the initial several hundred milliseconds after stimulation carry significant information about the target

To identify any possible trends between the Mutual Information Scores of 30 binned, discrete input features and the output target, we performed a line fitting to these scores. This was done to further understand any existing linear trends. The R² score for both Binary-Class and Multi-Class fitted lines were 0.85 and 0.87, respectively, and the slopes being -0.0047 and -0.0095, respectively. These negative slopes that are also evident in Fig. 8 indicates that the significant dependence of the target decreases with time post-stimulation.

Fig. 8

Mutual Information analysis for understanding feature significance in binary and multi-class classification. A Linear Trend analysis on Mutual Information Scores of 30 binned discrete features with the Binary-Class target, and (B) Linear Trend Analysis on Mutual Information Scores of 30 binned discrete features with the Multi-Class target. Each point indicates the mutual information score for a binned discrete input feature with target, and a fitted trend line showing the overall trend. The downward trends in both plots highlight that the significant dependence of the target on these features decreases over time, with the dependency diminishing as time approaches the end of the stimulation period (3000 ms). This analysis aids in understanding which features are most influential in predicting the target in both binary and multi-class scenarios

The feature importance analysis shown in Fig. 9 suggests that several input features are uniquely positioned to infer octopus arm movement, particularly within the first 100 ms period of the electrophysiology response.

Fig. 9

Decoding the top 10 Non-Linear dynamics via feature importance analysis. This analysis demonstrates the influence of different input features in model decision. It shows that for both (A) binary-class and (B) multi-class the type of stimulation and the timing of the feature is crucial

Computational models that predict movement

The dataset was segregated by utilizing the 80/20 Split method, where 80% (220 samples) of the data was used for cross-validating of the model, and 20% (55 samples) of the data was used for testing the finalized model. The Stratified k-Fold Cross-validation technique was applied with K of 10 folds.

The accuracies of 16 different classifying techniques are reported in Table 2. These average performance scores returned from our cross-validation tests. Additionally, despite the class imbalance in the 100 ms binary dataset, with class 1 consisting of 201 samples and class 0 consisting of 74 samples, the model still appears to be performant and able to handle the imbalance effectively based on the evaluation of F1 score results from cross-validation (Table 2), and test dataset (Table 3). The best model binary-class dataset was the Gradient Boosting Classifier, which achieved 88.64% accuracy. The best model for the multi-class dataset was the Extra Trees Classifier, which achieved an accuracy of 75.45%. These models also showed the highest F1 scores, describing the harmonic mean of precision and recall, which reflects the model's accuracy and consistency in identifying the correct category.

Table 2 Comparative analysis of machine learning classifier performance: cross-validation metrics for binary-class and multi-class tasks

Table 3 Evaluation of machine learning classifiers on test data: performance metrics across binary and multiclass classification

The performance metrics for the test data are presented in Table 3. The Gradient Boosting Classifier, identified as the best model through cross-validation on the Binary-Class dataset, was tested on 20% of previously unseen data and achieved an accuracy of 83.64%. This compares to the 88.64% accuracy obtained with the stratified k-Fold evaluation on the remaining 80% of the data. Similarly, the Extra Trees Classifier, which was the best model from cross-validation on a Multi-class dataset, was tested on 20% of unseen data and reached an accuracy of 72.73%, versus the 75.45% accuracy achieved with the stratified k-Fold evaluation on the remaining 80% of the data.

This analysis indicated that the Gradient Boosting Classifier could predict with 88.64% accuracy the movement in datasets composed of 30 bins of 100 ms each. The accuracy of movement prediction using 50 ms bins was also tested and compared to the 100 ms dataset. This dataset consisted of 60 continuous binned features, the categorical features, and the Target features. The results demonstrated that higher accuracy in movement prediction could be achieved with 100 ms compared to the 50 ms dataset, as shown in Table 4.

Table 4 Comparison of models across time frames

The confusion matrix (Fig. 10) is a method that allows computing a machine learning model accuracy, precision, recall, and overall ability to correctly classify instances, providing a detailed view of the types and frequencies of classification errors. These results suggest that the Gradient Boosting Classifier and the Extra Trees Classifier could predict type of movement with high accuracy.

Fig. 10

Confusion matrix analysis for binary and multi-class movement predictions. A Binary-class and (B) Multi-class confusion matrix analysis showing the correct and incorrect predictions of the type of movement based on the test dataset which is 20% of the entire data. The green and blue indicate the correct predictions. This analysis suggests the high accuracy of the models in predicting the type of movement

Computational models of movement to decode stimulation

We next considered a finer-grained analysis of the movements evoked from stimulation of the arm to determine relevant kinematic features beyond the 0, 1, and 2 movement class labels manually applied. Different stimulation locations elicited distinct behavioral responses in octopus arm movement. We examined the distribution of each kinematic metric previously defined, using both keypoint-derived features and proSVD identified bases. Figure 11 shows one example of how movement evoked from stimulation at the cord, versus at the distal or proximal regions (PD), has a significantly different distribution of the maximum angular velocity across each time period post-stimulation (see Appendix for complete table). Initial analyses considered the distal and proximal regions separately, but found no significant differences, and so all future analyses considered them as a single combined group.

Fig. 11

Different motion features evoked by stimulus location. A Histograms and overlaid kernel density estimation plots for the extracted maximum angular velocity in t0 for each video during (top) cord or (bottom) proximal or distal (PD) stimulation. B A principal component analysis was performed for all extracted behavior features. Selected features are shown here with higher loadings across the first three principal components (PC). The mean and maximum angular velocities and proSVD features across the first two time periods comprise much of the first two PCs, with the distal velocity in the last time period having the highest loading onto the third PC. C The distributions of this metric, maximum angular velocity, are significantly different across all post-stimulation time windows (for t0, t1, t2; p = 0.014, p = 0.001, p = 0.004, respectively, determined by a two-sample KS test)

Principal component analysis was conducted on all features across all time periods to identify which features might be playing the largest role (selected features shown in Fig. 11b). Immediately following stimulation, the angular speed of the distal part and the proSVD bases were the most significant kinematic features, with translation speed contributing more significantly in the later time periods. We found that the features with the highest loadings in PCA space also tended to produce significantly different feature distributions between stimulation types.

The above analyses all collapsed movement into three time periods, potentially missing relevant signals at finer time resolution. To characterize more complex behavior, and to establish real-time methods for future closed-loop work, we looked at what metrics and analyses we could do in the streaming setting, as fast as data could be collected. DLC live was able to generate keypoint inferences at rates of ~ 10–20 ms per image, or around 100 frames per second on average (Fig. 12a). proSVD could be run at extremely high frame rates, and may capture finer motion features than those that could be identified using keypoint tracking. We found that proSVD features showed some differences between mechanical and electrical stimulation types as a function of time post-stimulation (Fig. 12b), unlike earlier keypoint metrics. While less interpretable, we anticipate that including these types of unsupervised features alongside user-defined keypoints will produce the quantification needed to fully characterize the rich and complex behavior in the octopus repertoire.

Fig. 12

Real-time inference and time-varying features. a Inference speeds as a function of video resolution for both DeepLabCut live and proSVD during streaming analyses for different downsampled video resolutions. Error bars denote standard deviation (N = 10). b The L2-norm of the derivative of the first feature identified by proSVD plotted as a function of time post-stimulation, averaged across all trials. The low-dimensional representation of the evoked motion shows some different temporal features for each stimulus conditions

Discussion

The octopus’s extraordinary anatomy and physiology makes it especially attractive to uncover sensorimotor circuits that orchestrate behavior. The octopus’s nervous system is highly distributed, and much of the neural circuitry coordinating these behaviors is organized within the arms where it can monitor immediate, complex environmental feedback and adapt the arms’ movement accordingly (Sumbre et al. 2001, 2006, 2005; Levy et al. 2015). Identifying and understanding the neural signals that drive complicated motor output, such as those found in the octopus, are essential for the future development of rapidly-adapting and ultimately more human-like prosthetic arms (Sivitilli et al. 2022, 2023).

Here we describe high-temporal and spatial resolution results of single and multi-unit data, obtained from a detached, behaving octopus’s arm. The results show that the number of spikes occurring within the first 100 ms after stimulation can predict the movement response, whereas the greatest number of spikes were associated with lateral movement. Stimulation location was also a significant variable: The greatest number of spikes in the first 100 ms occurred in response to proximal stimulation, but distal stimulation evoked the greatest change in spike response over time. These results indicate that spike analysis can reveal fundamental principles of motor behavior. The statistical analysis of the stimulation location could predict with 58% accuracy what movement would be evoked. By adding more features, specifically the number of spikes as a function of time, to the ML algorithms, the accuracy of the prediction increased to 88.64%.

While the single and multiunit analysis by itself consisted of important information of behavior, we tested if using ML could further identify important features of octopus motor behavior. The models that were built on the ML results indicate that it was possible to predict whether an arm movement occurred with 88.64% confidence, and it was possible to predict with 75.45% confidence if this was a lateral arm movement or an arm movement with a curl. These levels of prediction values are in agreement with the range of reported accuracies in arm reaching predictions based on large-scale neural activity in monkey’s motor cortex (Tseng et al. 2019; Vaskov et al. 2018) and in humans (Batzianoulis et al. 2018). The accuracy of machine learning algorithms in predicting arm movements can vary depending on several factors, including the specific algorithm used, the quality and quantity of input data, the level of noise, and the complexity of the movement being predicted.

The findings presented here offer significant insights into the neural mechanisms underlying motor control, particularly in how distinct features of electrophysiological responses predict different types of movement. Our results demonstrate that the first 100 ms of the electrophysiological response is the most important predictor of movement type, which is consistent with previous studies suggesting that early neural activity plays a crucial role in motor planning and execution. Additionally, the stimulation location, particularly peripheral stimulation to the distal arm, emerged as a critical feature in our model. This aligns with existing evidence showing that distal stimuli elicit the most pronounced neural responses (Chang and Hale 2023), potentially due to an amplification mechanism in the distal regions of the arm that enhances sensitivity and reactivity. This finding opens the door for future research on the role of distal areas in fine motor control and the adaptive behaviors they facilitate, such as searching and exploring.

A key discovery from our motion analysis is that the same stimulation does not always lead to the same movement, highlighting the complexity of motor control that relies on both central and peripheral neural circuits (Zullo et al. 2009, 2019; Hochner et al. 2023). The probability distribution of movements, rather than a singular response, further supports the idea that motor control is governed by a dynamic and adaptive network. This insight points to the presence of decentralized intelligence in the motor system, where timing, location, and environmental factors collectively influence motor circuit activation. Such flexibility allows organisms like the octopus to exhibit a high degree of adaptive behavior beyond simple reflexive responses.

Traditionally, it has been believed that all motor computation occurs in the brain and central nervous system. However, our findings challenge this view, providing evidence that a significant amount of computation can occur outside the brain—specifically in the peripheral nervous system and within the arm itself. This peripheral computation, in conjunction with central processes, contributes to the generation of meaningful and complex movements. This discovery highlights the role of the arm and distal regions in motor control, suggesting that these areas play a more active role than previously thought in shaping movement patterns. This raises important questions about how we model motor control and the potential for decentralized systems to influence behavior in ways that were once considered solely the domain of the central nervous system.

These findings also highlight a fundamental challenge in current models of motor control: there may not be a single, unified model capable of predicting all aspects of movement in the context of intentional motion. Intentional movements likely require computations that involve multiple brain regions and feedback loops that are yet to be fully identified and explored. As our results show, the complexity of movement is influenced by many factors, including timing, location of stimulation, and environmental context. These factors likely engage distributed neural networks and dynamic feedback systems that evolve over time. The decentralized nature of motor control complicates the creation of a one-size-fits-all model, emphasizing the need for further research to map out the interconnected neural pathways and feedback mechanisms that contribute to the execution of intentional movement.

We have recently used a set of five reflective markers that were adhered to the octopus skin, in order to describe and quantify the overall posture of an awake, swimming octopus (Weidig et al. 2024). Three postures that were defined as straight, simple bending, and complex bending, and were analyzed in 3D using curvature and plane orientation methods. The results showed that this novel kinematics approach was successful in understanding octopus posture. However, this approach is limited by the number of markers that could be attached to the octopus’s arm and the physical constrains they might induce.

Nevertheless, state-of-the-art computer vision and machine learning tools could provide quantification of kinematic features based on video recordings alone (Mathis and Mathis 2020; Pereira et al. 2020; Goodwin et al. 2024). Deep learning and other machine learning methods could also learn features that human eyes do not see, but may be significantly correlated with neural firing or stimulation patterns (Calhoun et al. 2019; Schneider et al. 2023; Syeda et al. 2024). Employing transfer learning, deep neural networks, and dimension reduction as described here, we aimed to gain real-time insights into octopus arm movements and how their motor circuitry produces rich movement types.

There are limitations in current methods for accurately locating points of interest. In the octopus arm videos considered here, several display minimal to no movement except for mechanical adjustments made by the experimenter. The predominance of single-instance movement in most videos, with shorter clips, effectively limits our dataset for training the DLC model on complicated movements. Additionally, the non-planar motion of the arm at times poses a challenge for accurate tracking, requiring new tracking strategies to capture the full complexity of octopus arm movements.

Our selection of kinematic parameters was inspired by a study on locomotion using zebrafish larvae (Marques et al. 2018), but their effectiveness for the octopus arm is not straightforward due to distinct ethology and experimental conditions. Unlike zebrafish, the octopus arm lacks a zero-angle "tail" at rest. Zebrafish data were made positionally consistent through affine transformations and background removal, a step not directly applicable to octopus arm data due to sample size limitations, hindering common clustering techniques. Assessing the effectiveness of unsupervised clustering to identify key features also proves challenging in the absence of ground truth labels to gauge cluster accuracy and precision. Overall, improvements could be achieved through higher resolution videos, camera stability without flickering, incorporation of multiple stable camera angles, precise manual annotation, and a larger volume of data.

In conclusion, this study provides valuable insights into the complexities of motor control, particularly in the context of octopus arm movements. Through meticulous video analysis, electrical and mechanical stimulation, and kinematic feature extraction, we have gained a deeper understanding of the key factors influencing arm motion, such as angles, angular speed, and absolute speed. These findings highlight the intricate interactions between neural circuits and peripheral feedback mechanisms, suggesting a decentralized model of motor control that goes beyond traditional views. Machine learning techniques offer great potential in unraveling these complexities by simulating how various neural systems interact and predict movement patterns. By incorporating critical features, such as stimulation timing and location, into machine learning models, we can develop adaptive algorithms that simulate the underlying motor processes and enhance the precision of brain-machine interfaces, prosthetics, and neural stimulation therapies. Ultimately, this research not only advances our understanding of motor behavior and neural circuitry but also lays the groundwork for designing systems that replicate natural, adaptive movement. These advancements are crucial for the development of more sophisticated, flexible control systems in both biological and artificial contexts, moving us closer to creating technologies that mimic the precision and adaptability of natural motor control.

Appendix

Table 5 K-S Significance Test on PD (Proximal–Distal v/s Cord)