A spiking neural network model for fractional proprioceptive encoding of limb posture and movement in insects

TL;DR

This paper introduces a spiking neural network model that explains how insects encode limb posture and movement using fractional proprioceptive signals.

Contribution

A novel hierarchical spiking neural network model for distributed computation of posture and movement in insects.

Findings

Adaptive Exponential Integrate-and-Fire neurons model phasic-tonic encoding of joint angles by proprioceptive afferents.

Interneurons accurately encode joint angle and velocity across the full working range.

The hierarchical model captures complex movement patterns from single joints to whole-body posture.

Abstract

Proprioception is key to all behaviours that involve the control of force, posture or movement. Computationally, many proprioceptive afferents share three features: First, their strictly local encoding of stimulus magnitudes causes range fractionation in sensory arrays. As a result, encoding of large joint angle ranges requires convergence of afferent information onto first-order interneurons. Second, their phasic-tonic response properties lead to fractional encoding of the fundamental sensory magnitude and its derivatives (e.g., joint angle and angular velocity). Third, the distribution of disjunct sensory arrays across the body implies that complex movements involve information from multiple joints or limbs. The present study proposes a multi-layer spiking neural network for distributed computation of whole-body posture and movement. The first part of the study models strictly local, phasic-tonic encoding of joint angle by proprioceptive hair field afferents by use of Adaptive Exponential Integrate-and-Fire neurons. Fractionally encoded afferent information about single-joint posture and movement converges on two types of first-order interneurons, tuned to encode either joint angle or velocity across the entire working range with high accuracy. In velocity-encoding interneurons, spike rate increases linearly with angular velocity. The companion paper exploits this distributed position/velocity encoding in second- and third-order interneurons, using combinations of two or three position/velocity inputs from disjunct arrays. The encoding properties of all interneuron layers are evaluated with experimental data on whole-body kinematics of unrestrained stick insect locomotion, comprising concurrent joint angle time courses of \documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$6\times 3$$\end{document} leg joints. The hierarchical model allows increasingly complex encoding of posture and movement, from angular velocity of a single joint, to movement cycle phases of an entire limb, to parameters of overall body posture.

Figures

Distributed proprioceptive encoding of leg posture and movement A. The six legs of the Indian stick insect *Carausius morosus* have similar size and structure. The dataset used here includes joint angles of the thorax-coxa, coxa-trochanter and femur-tibia joints, referred to as the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}$ , $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepack

Bi-directional hair field and extended bi-directional hair field. A. The bi-directional hair field models pairs of hair rows on opposite sides of the joint and working-ranges below (red) and above (blue) the neutral angle $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${90}^{\circ }$$\end{document}$ . In the idealised hair row, hair angle is a function of joint angle with parameters $\documentclass[12pt]{minimal}

Phasic-tonic spike response of the real and modeled mechanosensory neuron. A. Spike response of mechanosensory afferents from the lateral scapal hair plate in the American cockroach (Okada and Toh [2001](#CR47)). The plot shows the spike rate over time for a tactile hair deflected by a ramp-and-hold function with a constant velocity of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${60.1} \frac{\circ }{{\textrm{s

Position IN results. A. Raster plot displaying the response of a hair plate ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h} = 50$$\end{document}$ ) to a 5-second movement sequence of the $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \be

Velocity IN results. A. For velocity INs, the output spike rate increases linearly with angular velocity. $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document}$ and $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\t

Full text

Introduction

Proprioception is the ability to perceive body posture, movement, and load, a ubiquitous ‘sixth sense’ in all mobile organisms (Tuthill and Azim 2018). It relies on information provided by mechanosensory neurons, known as proprioceptors. In insects and other arthropods, proprioceptors are distributed throughout the entire musculo-skeletal system, most of them embedded in or attached to the cuticle (McIver 1985). Multiple types of mechanoreceptors encode physical magnitudes (Tuthill and Wilson 2016), as different as force or strain (Campaniform sensilla: e.g., Pringle 1938; Chapman 1965; Hofmann and Bässler 1986; reviewed in Zill et al. 2004), joint angle or angular velocity (Chordotonal organs: e.g., Hofmann et al. 1985; Mamiya et al. 2018; reviewed by Field and Matheson 1998; hair fields: e.g., Pringle 1938; Wong and Pearson 1976; Schmitz 1986) with a high degree of convergence of afferent information from distinct sensors (Stein and Schmitz 1999; Schmitz and Stein 2000; Gebehart et al. 2022). Quite generally, all proprioceptor types are distributed across the body, with disjunct sensory arrays on multiple limbs and limb segments (e.g. Markl 1962) with consistent structure-function relationships across different species (e.g. Virdi and Sane 2024) and across limb types as different as walking legs and antennae (Krishnan and Sane 2015).

Here, we propose a computational model for hierarchical processing of distributed proprioceptive information about limb and body kinematics. As key variables of internal representations of posture and movement, we focus on joint position and velocity. In contrast to many models available (see below), we use a spiking neural network (SNN), to reflect the spiking nature of all prorioceptive afferents (Mill 1976) and many of their downstream interneurons (Burrows 1996). An important feature of SNNs is their similarity to networks of the central nervous system (CNS), integrating neural and synaptic states, temporal dynamics, and the generation of action potentials (Yamazaki et al. 2022).

The main objective of our study is to encode increasingly complex information about posture and movement, spanning the range from single joint angles, to the movement of an entire limb and, finally, to whole-body posture. We do so using a feed-forward network, thus neglecting the involvement of proprioceptive information in dynamic feed-back control, and focusing on internal representation and state estimation instead (Dallmann et al. 2021).

This first part of this study addresses two major computational challenges of peripheral proprioceptive encoding in general: (i) fractional encoding of a stimulus variable and its derivatives, and (ii) range fractionation in sensory arrays.

The challenge of fractional encoding is a computational consequence of spike rate adaptation, as evident in the phasic-tonic spike frequency response to step changes in the stimulus magnitude. For example, chordotonal organs and hair fields from different species and different limbs display a mix of tonic (slowly adapting) and phasic (rapidly adapting) response components to step changes in joint angle or angular velocity (Zill 1985; Hofmann et al. 1985; Newland et al. 1995; Okada and Toh 2001). This is thought to be caused by both the mechanics of the surrounding or embedding structure (e.g., Barth 2019), and by the properties of the sensory neuron (e.g., French et al. 2002). As a consequence, the afferent spike rate encodes the input magnitude itself (e.g., joint angle) along with its derivatives (e.g., angular velocity). Models of specific sensorimotor pathways have exploited this fractional encoding for computational purposes in individual interneurons (INs) (Jones and Gabbiani 2012) or for PD control of steering (Cowan et al. 2006).

Apart from models concerning the cellular mechanisms underlying spike rate adaptation (e.g., French et al. 2002), computational models of phasic-tonic proprioceptor responses have been based on different levels of description. For example, Cocatre-Zilgien and Delcomyn (1999) modeled the afferent spike rate of campaniform sensilla by means of a two-stage stimulus-response function, where the first stage captured the tonic component as a hyperbolic function of strain (in analogy to vertebrate mechanoreceptors: Loewenstein 1961), and the second stage implemented phasic adaptation by means of a power law, as previously proposed for mechanoreceptor adaptation in cockroaches (Chapman and Smith 1963; French 1984). Applying a similar approach to hair fields, “total afferent activity” of the stick insect trochanteral hair field has been simulated using a phasic-tonic function of spike rate on joint angle (Dean 1985). In contrast, sensory array models involve multiple parallel receptor models. For example, Ache and Dürr (2015) applied cascaded low-pass and high-pass filter blocks to model different stages of proprioreceptive encoding of antennal position and velocity in stick insects. Similarly, Szczecinski et al. (2021) modeled phasic-tonic changes in a strain-sensitive campaniform sensillum. While this kind of analog stimulus-response functions may be combined with stochastic spike generators to generate time sequences of spike time events (e.g., Gollin and Dürr 2018), there is a lack of a spiking proprioceptor model that generates spike trains directly through subthreshold membrane potential dynamics. Our present model does exactly this, employing the mathematically compact Adaptive Exponential Integrate-and-Fire (AdEx) (Brette and Gerstner 2005) for rate adaptation of a spiking neuron. Moreover, we exploit the fractional encoding properties of this AdEx model to show that downstream INs can decode both fractionally encoded magnitudes - in our case joint angle and joint angle velocity - with high precision. As a particular example, we attempt to model linear velocity encoding across a wide range of input velocities, as described for first-order INs of the antennal mechanosensory pathway (Ache et al. 2015).

The second computational challenge concerns the fact that, quite generally, mechanotransduction encodes local forces and/or deformations. Therefore, a common feature of all proprioceptive afferents is their strictly local encoding of sensory magnitudes. As a consequence, proprioceptive organs are sensory arrays that show range fractionation of their overall receptive field (Matheson 1992). Computationally, this requires convergence of multiple afferents to obtain a first-order IN with a suitably large receptive field. In a recent biorobotics approach, Zadokha and Szczecinski (2024) proposed range-fractionated sensory arrays to infer foot position from a set of joint angles. In their approach, the sensory input was provided by a juxtaposed series of analog receptors with Gaussian receptive fields. The synaptic weights between these and 1st- and 2nd-order INs were tuned to encode an analog quantity of choice, for example the robot’s foot height above ground. In contrast to these static analog joint angle receptors, our study proposes an array of direction-selective (non-linear) spiking afferents that generate time-varying, phasic-tonic spike responses, with fractional encoding of both joint angle and angular velocity. To do so, our 1st-order IN model uses a simple Leaky Integrate-and-Fire (LIF) neuron model along with appropriate spatial convergence of the sensory input array.

For quantitative comparison with published experimental data on posture and movement of insect limbs, we chose to model the phasic-tonic spike response of hair field afferents of a cockroach, while evaluating downstream position and velocity encoding against whole-body kinematics data from walking stick insects. With regard to our objective and the two computational challenges explained above, we use detailed electrophysiological data on phasic-tonic single-hair spike trains in response to controlled changes of position and velocity at the antennal scapal hair plate of the American cockroach Periplaneta americana (Okada and Toh 2001). Downstream encoding of position and velocity of natural limb movements was evaluated with motion capture data on unrestrained walking and climbing stick insects Carausius morosus (Theunissen and Dürr 2013). Stick insects bear proprioceptive hair fields at all limb bases (antennae: Krause et al. 2013; legs: Wendler 1964; Cruse 1976; Fig. 1A) where they are located on the cuticle near joint membranes (Fig. 1B). A detailed description of size, number and arrangement of hair fields at the stick insect thorax-coxa and coxa-trochanter joints can be found in Schmitz (1985b). Scanning electron micrographs thereof may also be found in Cruse et al. (2009). When deflected, the hair acts as a lever arm, exerting a force on a sensory neuron dendrite (Tuthill and Wilson 2016), typically causing a phasic-tonic response (Newland et al. 1995; Okada and Toh 2001). In stick insects, ablation of individual hair fields affects leg positioning (Wendler 1964; Cruse et al. 1984) and joint angle control (Kemmerling and Varju 1982; Schmitz 1986) in standing and walking animals (Schmitz 1985a), affecting swing height during walking (Theunissen et al. 2014) and antennal inter-joint coordination (Krause et al. 2013) alike.

We use a hierarchichal, four-layered SNN to accurately encode joint, limb and whole-body kinematics from a set of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$18\times 2$$\end{document} proprioceptive sensory arrays with phasic-tonic response properties (Fig. 1D), as evaluated by experimental data comprising concurrent joint angle time courses from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$6\times 3$$\end{document} joint angles (Fig. 1C). This first of two companion papers addresses the fundamental problem of fractional encoding of joint angle and angular velocity in the presence of range fractionation in disjunct sensory arrays (SNN layers 1 and 2). The second part then builds on the distributed internal representation of limb kinematics to test the decodability of higher-order parameters like intra-leg movement primitives as characteristic signatures of particular step cycle phases (SNN layer 3), and of the body pitch angle from spatial inter-leg coordination (SNN layer 4; van der Veen et al. 2026).

The remainder of this work is structured as follows: The methods section provides an introduction to the dataset and a detailed explanation of the network architecture, describing the methodology layer by layer. The result section presents our findings for each layer of the network. The discussion section provides a comparative analysis with the relevant scientific literature, explores the strengths and weaknesses of the proposed model, and ends with an outlook on future research.

Fig. 1. Distributed proprioceptive encoding of leg posture and movement A. The six legs of the Indian stick insect Carausius morosus have similar size and structure. The dataset used here includes joint angles of the thorax-coxa, coxa-trochanter and femur-tibia joints, referred to as the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} joints, respectively. Red arrows indicate locations of proprioceptive hair fields at the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} joints of the legs (labeled Cx and Tr, respectively) as well as the head-scape and scape-pedicel joints on the antenna (labeled Sc and Pd, respectively). B. Schematic representation of a proprioceptive hair field on the cuticle near a joint, before and during the joint membrane folds over the hair field during a change in joint angle. A change in joint angle causes hair deflection. As the joint angle increases, the number of hairs deflected and the deflection angle per hair also increase. Each hair acts as a lever arm connected to a dendrite of a sensory neuron. Mechanotransduction channels open, allowing the sensory neuron to spike in proportion to the hair deflection (Tuthill and Wilson 2016). C. Example time course of the joint angle for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} joint of the right anterior leg (R1). D. Proposed SNN architecture for distributed proprioception of limb kinematics and body posture. The joint angle time courses \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} for each of the six legs are converted into hair angle time courses (not shown) that are sensed by one mechanoreceptor per hair in Layer 1. Solid and open circles refer to antagonistic arrangement of two hair fields per joint. In Layer 2, two velocity (vel +/-) and position (pos +/-) INs per joint encode the posture and movement. In the companion paper, Layer 3 integrates converging information from all joints per limb, and Layer 4 integrates converging information from all limbs

Methods

Dataset

The experimental data used in this study was originally collected to study whole-body kinematics of walking and climbing stick insects, comparing related species with different body morphology (Theunissen et al. 2015) and characterizing distinct step classes (Theunissen and Dürr 2013). Stick insects of the species Carausius morosus (de Sinéty, 1901) have six legs with fairly similar morphological structure (Fig. 1A). Each leg comprises a short basal coxa, a fused trochantero-femur, a long and thin tibia, and distal tarsus with five tarsomeres. The dataset contains motion capture data on three joint angles per leg. Two of these joint angles correspond to joints monitored by proprioceptive hair fields: Hair fields on the coxa (red arrows labelled Cx) measure protraction-retraction movements of the thorax-coxa joint, denoted by a blue \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} in Fig. 1A (Wendler 1964; Cruse 1976); The trochanteral hair field (red arrows labelled Tr) measures levation-depression movements of the coxa-trochanter joint, denoted by a green \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} in Fig. 1A (Theunissen et al. 2014). The dataset also contains extension-flexion movements about the femur-tibia joint, denoted by a red \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} in Fig. 1A. Throughout this work, these three joints are referred to as the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} joints, respectively. The front, middle, and hind legs are labeled as 1, 2, and 3 for the right (R) and left (L) sides. Fig. 1A also indicates the location of antennal hair fields at the head-scape and scape-pedicel joints on the antenna (labelled by Sc and Pd for scape and pedicel, respectively). Since the dataset does not comprise antennal joint angles, the current study focuses on leg proprioception only. However, given similar function of all proprioceptive hair fields, the spiking network proposed in this work could potentially be expanded to include antennal kinematics. Generally, we assume that the joint angle encoding mechanism follows the same principle for all mentioned joints and all limbs.

The dataset comprises complete body kinematics of unrestrained walking and climbing stick insects. Nine specimens walked freely on a horizontal walkway measuring 40 mm in width and 490 mm in length. In this work, we focus on trials where the animals encountered a flat surface, whereas the companion paper expands to include climbing trials with two stairs (Theunissen and Dürr 2013). A marker-based motion capture system (Vicon MX10) was employed, using eight infrared cameras capturing 200 frames per second, to track markers attached to the head, thorax, and all six legs of the insect. The captured marker trajectories were used to reconstruct the time courses of the three mentioned leg joint angles \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} .

Spiking neural network (SNN) architecture

A schematic of the proposed SNN architecture is shown in Fig. 1D. The temporal evolution of each experimentally obtained joint angle is converted into a set of hair deflection angles, corresponding to the number of hairs in the hair field. Each hair has a unique receptive field that is arranged in sequence with receptive fields of adjacent hairs. The combination of all receptive fields results in sensitivity for the entire working range of the joint. In analogy to real proprioceptive hair fields (Fig. 1B), each hair deflection is converted into an electric current that is integrated by a single mechanosensory neuron per hair. The mechanosensory neuron consists of an AdEx model (Brette and Gerstner 2005) and provides a nonlinear phasic-tonic response to a step increase in stimulus magnitude (here: joint angle). This is due to an adaptation mechanism that modulates the relationship between the input current and the output spike rate. The resulting spike trains of all sensory neurons per hair field converge on a set of four first-order INs per joint: two position and two velocity INs. One position INs encodes negative joint displacement relative to rest (denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pos_-$$\end{document} ), while the other encodes positive joint deflection relative to rest ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pos_+$$\end{document} ). Similarly, the velocity INs encode a change in joint angle into either the positive (low joint angle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow$$\end{document} high joint angle) or negative (high joint angle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow$$\end{document} low joint angle) direction, denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_+$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_-$$\end{document} , respectively. For simplicity, each leg of our model has three joints, each with its own hair plate containing two hair fields with the same proprioceptive encoding mechanism, assuming that key features such as receptive field size, linear range fractionation with equal sensitivity, and phasic-tonic response time courses are the same at all 18 leg joints. Accordingly, the model uses 3 sets of two hair field implementations per leg, totalling \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3 \times 2 \times N_\text {h}$$\end{document} mechanosensory neurons that converge onto \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3 \times 2 = 6$$\end{document} position INs and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3 \times 2 = 6$$\end{document} velocity INs. Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h}$$\end{document} represents the number of hairs in a hair field. Both IN types are modelled by a simplified LIF model, providing a linear relation between an input current and output spike rate.

The proposed architecture was implemented in Python version 3.9. The simulations were conducted on a system equipped with 16 GB of RAM and an AMD Ryzen 5600x processor. The differential equations that describe the behaviour of spiking neurons and synapses were solved over time using the backward difference method and a time step of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$dt = {0.25} \hbox { ms}$$\end{document} . The Zenodo repository containing Python-Jupyter notebooks can be accessed via the following link.

Layer one: Hair plate

Hair field

It is important to note that while proprioceptive hair fields are present at the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} joints of each stick insect leg, they are absent at the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} joint. Despite this fact, hair field models were applied to all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} joints as well because mechanosensory neurons of the femoral chordotonal organ, i.e., the sensory organ encoding the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} joint angle, share important encoding features with mechanosensory neurons of proprioceptive hair fields, such as sensitivity to joint angle, joint angle velocity (Hofmann et al. 1985), and range fractionation (Matheson 1992; Ache and Dürr 2013).

In each hair field, the time course of the respective joint angle is transformed into a set of hair deflection angles, corresponding to the number of hairs on a hair field, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h}$$\end{document} . Each hair features a distinctive receptive field arranged in series with its neighboring hairs. The receptive field of each hair is defined as the range of joint angles within which the sensillum is sensitive to changes in deflection. For simplicity, the collective contribution of all hairs spans the entire possible range of the joint angle. Moreover, the receptive field size and the spacing between hairs are uniform for all hairs within a given hair field, deviating from the variation observed in biological hair fields (Pringle 1938). It is also assumed that there is a degree of overlap between receptive fields. This accounts for the fact that many hair fields are not arranged in regular hair rows but form patches (or plates) of hairs with overlapping receptive fields. Generally, it is unlikely that full deflection of one hair coincides perfectly with the onset of deflection of the next hair. Finally, it is assumed that hair deflection is linearly proportional to the joint angle and the range is bound to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[{0}^{\circ }, {90}^{\circ }]$$\end{document} . If the joint angle falls below or exceeds the proprioceptor receptive field, the hair angle is either not deflected at all ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${0}^{\circ }$$\end{document} ) or fully deflected ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${90}^{\circ }$$\end{document} ), respectively. This results in the following relation between joint angle ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta$$\end{document} ) and hair angles ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi _{ij}$$\end{document} ) for hair i and hair field j:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \phi _{ij}(\theta )= {\left\{ \begin{array}{ll} 0& \text {if } \theta < \theta ^\text {(rf0)}_{ij}\\ 90& \text {if } \theta> \theta ^\text {(rf90)}_{ij}\\ \frac{90(\theta -\theta ^\text {(rf0)}_{ij})}{\theta ^\text {(rf90)}_{ij} -\theta ^\text {(rf0)}_{ij}}& \text {otherwise} \end{array}\right. } \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^\text {(rf0)}_{ij}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^\text {(rf90)}_{ij}$$\end{document} are the lower and upper receptive field edges, respectively, defined as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \theta ^\text {(rf0)}_{ij}&= \frac{\theta _{j}^{(\text {max})}-\theta ^{\text {(min)}}_{j}}{N_\text {h}}(i-1) + \theta ^{\text {(min)}}_{j} - \frac{\theta ^{(\text {ol})}_j}{2} \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \theta ^\text {(rf90)}_{ij}&= \frac{\theta _{j}^{(\text {max})}-\theta ^{\text {(min)}}_{j}}{N_\text {h}}i + \theta ^{\text {(min)}}_{j} + \frac{\theta ^{(\text {ol})}_j}{2} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{j}^{(\text {max})}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^{\text {(min)}}_{j}$$\end{document} represent the maximum and minimum of the the joint angle sensitivity for hair row j, respectively. In the proposed network, these parameters are set as the maximum and minimum joint angles attained by the corresponding joint. Therefore they are not equal for all hair fields. The receptive fields of the outer hairs are manually set to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{\text {rf0}}^{(1j)} = \theta ^{\text {(min)}}_{j}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _{\text {rf90}}^{(Nj)} = \theta _{j}^{(\text {max})}$$\end{document} . The overlap between two adjacent receptive fields is defined by the parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^{(\text {ol})}_j$$\end{document} and is equal for all hair fields. Similar to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h}$$\end{document} , which is also consistent across all hair fields.

Bi-directional hair fields

Hair fields such as the antennal (Krause et al. 2013) or coxal (Wendler 1964) hair fields of stick insects are often arranged in opposing pairs. Accordingly, our model divides proprioceptive hair fields at each joint into two complementary and opposing subgroups. The parts of these bi-directional hair fields (Fig. 2) have their working-range either in the upper or lower half of the joint angle working-range. The halfway point of the joint working-range is called the ‘resting angle’ or ‘neutral angle’, at which no hair is fully deflected. The two hair fields deflect due to increasing and decreasing joint angle relative to the resting angle, respectively. Hair fields with opposing deflection sensitivity require a modification in Eq. (1):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \phi _{ij}(\theta )= {\left\{ \begin{array}{ll} 90& \text {if } \theta < \theta ^\text {(rf0)}_{ij}\\ 0& \text {if } \theta> \theta ^\text {(rf90)}_{ij}\\ 90\left( 1-\frac{\theta -\theta ^\text {(rf0)}_{ij}}{\theta ^\text {(rf90)}_{ij} -\theta ^\text {(rf0)}_{ij}}\right) & \text {otherwise} \end{array}\right. } \end{aligned}$$\end{document}

Moreover, the calculations for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^\text {(rf0)}_{ij}$$\end{document} (Eq. (2)) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^\text {(rf90)}_{ij}$$\end{document} (Eq. (3)) are interchanged.

Fig. 2A illustrates multiple hair angles in relation to the joint angle for a hypothetical scenario. The calculations for positively and negatively oriented hairs are determined using Eqs. (1) and (4), respectively. The positively oriented hair field is sensitive from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${88}^{\circ }$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${180}^{\circ }$$\end{document} while the negatively oriented hair field is sensitive from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${92}^{\circ }$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${0}^{\circ }$$\end{document} . Incorporating an overlap into the hair plate was intended to maintain low but non-zero spike rates at the resting position for the hair field. The overlap between the hair fields ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^{\text {(olhf)}}$$\end{document} ) at the resting position ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${90}^{\circ }$$\end{document} in the hypothetical scenario) is defined as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \theta ^{\text {(olhf)}} = \theta _1^{(\text {max})} - \theta _2^{(\text {min})} = {92}^{\circ } - {88}^{\circ } = {4}^{\circ } \end{aligned}$$\end{document}

In the proposed network, the parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^{\text {(olhf)}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^{(\text {ol})}_j$$\end{document} are set to the same value.

Fig. 2. Bi-directional hair field and extended bi-directional hair field. A. The bi-directional hair field models pairs of hair rows on opposite sides of the joint and working-ranges below (red) and above (blue) the neutral angle \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${90}^{\circ }$$\end{document} . In the idealised hair row, hair angle is a function of joint angle with parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _1^{(\text {min})} = {0}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _1^{(\text {max})} = {92}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^{(\text {ol})}_1 = {4}^{\circ }$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h} = 5$$\end{document} (solid red lines), and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _2^{(\text {min})} = {88}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _2^{(\text {max})} = {180}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^{(\text {ol})}_2 = {4}^{\circ }$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{\text {h}} = 5$$\end{document} (dotted blue lines). The positively and negatively oriented hair angles are calculated using Eq. (1) and Eq. (4), respectively. The resting angle is at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${90}^{\circ }$$\end{document} . B. Extended bi-directional hair plate: This hypothetical scenario is modified from panel A as follows: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{\text {h}} = 10$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _1^{(\text {max})} = {180}^{\circ }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _2^{(\text {min})} = {0}^{\circ }$$\end{document} . Overlap ranges have been omitted for clarity.

Adaptive Exponential Integrate-and-Fire (AdEx) model

A deflected hair acts as a lever arm and applies force to the tip of the sensory neuron’s dendrites, opening mechanotransduction channels and generating a current (Thurm 1965). To reflect this in the model, the hair angles calculated in Eqs. (1) and (4) were multiplied by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}/{\phi } =$$\end{document} 10–150 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{{pA}}{\circ }$$\end{document} , yielding currents in the nA range. This current flows into a single sensory neuron, converting the input current into a spike train (Tuthill and Wilson 2016). To model the mechanosensory neuron dynamics of insect hair fields, we used empirical data on the dynamics of antennal proprioceptive hair field afferents characterized by Okada and Toh (2001) , illustrated in Figs. 3A,B. These experiments were carried out on the antenna of the American cockroach Periplaneta americana. Therefore, the mechanosensory mechanisms can be assumed to be similar between species and location, although particular parameters may need to be adapted to reflect differences in biology. Multiple neuron models were found suitable to replicate these non-linear spiking dynamics (Izhikevich 2004). While carefully considering biological accuracy, implementation costs, and potential spiking dynamics, the Adaptive Exponential Integrate-and-Fire (AdEx) model was chosen as the mechanosensory neuron model (Brette and Gerstner 2005). The AdEx model dynamics are governed by two ordinary differential equations (ODEs):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&C\frac{dV(t)}{dt} = I(t)-g_\text {L}(V(t)-E_\text {L}) + \nonumber \\ \quad&+ g_\text {L}\Delta _\text {T} \exp \left( \frac{V(t)-V_\text {T}}{\Delta _\text {T}}\right) -w(t), \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&\tau _\omega \frac{dw(t)}{dt} = a(V(t)-E_\text {L})-w(t), \end{aligned}$$\end{document}

where V(t) represents the membrane potential, C denotes the capacitance, I the input current, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_\text {L}$$\end{document} the leak conductance, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_\text {L}$$\end{document} the leak reversal potential, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta _\text {T}$$\end{document} the slope factor, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_\text {T}$$\end{document} the threshold voltage, w(t) the adaptation variable, a the adaptation coupling factor, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _w$$\end{document} the adaptation time constant (Brette and Gerstner 2005).

In a real neuron, an action potential (spike event) occurs due to depolarization. In the AdEx model, when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V> V_\text {t}$$\end{document} , a spike is recorded, and the timestep is noted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^f$$\end{document} . At this timestep, V(t) is reset to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_\text {L}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega (t)$$\end{document} increases by a constant b.

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {if} \hspace{2mm} t = t^f \hspace{2mm}&\text {then} \hspace{2mm} V(t^f) \rightarrow E_\text {L} \nonumber \\&\text {and} \hspace{2mm} \omega (t^f) \rightarrow \omega (t^f) + b \end{aligned}$$\end{document}

The second term in Eq. (6) on the right-hand side (RHS) represents the passive membrane function, implemented as a leakage mechanism that allows the membrane voltage to return to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_\text {L}$$\end{document} , in the absence of an applied current. In a biological neuron, this leakage term represents the random diffusion of ions across the membrane. The third term allows for the active (or voltage-dependent) membrane properties, whereby the membrane voltage triggers an exponential spike if it exceeds the threshold voltage, resulting in a transient, positive overshoot of the membrane potential. The fourth term on the RHS corresponds to the adaptation term, which modulates the active membrane properties by means of a current flowing out of the membrane. Consequently, an increase in w(t) results in a decrease in the sensitivity of V to a stimulus. The adaptation variable increases when the membrane voltage exceeds its resting state (Eq. (7)) or after a spike event (Eq. (8)). The second term on the RHS of Eq. (7) allows the adaptation variable to converge to zero in the absence of spikes, thereby resetting the adaptation process. Therefore, adaptation reduces the neuron’s sensitivity to a stimulus due to a prolonged or repeated stimulus. For reference, the response of the AdEx model to a constant current stimulus is shown in supplementary Fig. S1.

Layer two: Position interneurons (INs)

The overall activity of the hair plate could be decoded as the temporal evolution of joint angles by an IN that integrates spikes from all proprioceptors in a hair field. Due to the binary nature of the hair field (Fig. 2A), two position neurons were linked to the sensory neurons of negatively and positively oriented hair fields, respectively. With this arrangement, only one position neuron becomes active depending on whether the joint is within the negative or positive working range relative to its resting position, similar to the position encoding INs identified by Ache and Dürr (2013) .

Leaky Integrate-and-Fire (LIF) model

Due to its integration capabilities, simplicity, and small parameter set, the LIF model was chosen as the position IN model (Izhikevich 2004). In this model, a pre-synaptic spike that occurs at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^\text {pre}$$\end{document} increases V(t) by a synaptic weight \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document}  (Burkitt 2006):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} if \hspace{2mm}t=t^{pre}\hspace{2mm} then \hspace{2mm} V = V + \omega \end{aligned}$$\end{document}

The dynamics of the LIF are as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{dV(t)}{dt} = -\frac{(V(t)-E_\text {L})}{\tau } \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} is the decay time constant. If the voltage V(t) exceeds the threshold voltage \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_\text {T}$$\end{document} , a spike time was recorded as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^f$$\end{document} and the voltage was reset to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_\text {L}$$\end{document} at this timestep:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} if \hspace{2mm}t=t^f\hspace{2mm} then \hspace{2mm} V(t^f) \rightarrow E_\text {L} \end{aligned}$$\end{document}

The only tunable parameters are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} . For reference, the response of the LIF model to a constant frequency spike train stimulus is provided in supplementary Fig. S2.

Layer two: Velocity interneurons (INs)

First-order mechanosensory INs can also encode changes in joint angle. For example, the spike rate of a movement-sensitive IN of the antennal mechanosensory system increases linearly with joint velocity but remains inactive when the joint is stationary (Ache et al. 2015). To replicate this kind of velocity encoding, velocity INs are modeled to spike whenever the joint angle transitions from the receptive field of one hair to the next. To accomplish this, we could leverage the phasic response of the sensory neurons.

Modified hair field distribution

Similar to the position INs, the movement layer features two velocity INs per joint. The INs fire during an increase ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_+$$\end{document} ) or decrease ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_-$$\end{document} ) of joint angle, respectively. However, to encode velocity across the entire working range, both neurons need to have sensitivity throughout the complete joint range of motion. Owing to the directional selectivity of each hair’s response, this is not achievable with the opposing arrangement of the simple bi-directional hair field. Consequently, the scenario shown in Fig. 2A was extended by including 5 additional hairs and expanding the joint range for both hair fields: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h} = 10$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _1^{(\text {max})} = {180}^{\circ }$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _2^{(\text {min})} = {0}^{\circ }$$\end{document} . This modification is hypothetical in that it could either account for a response to ‘release of deflection’, or for a fraction of hairs being fully deflected at the rest position of the joint. Both of these properties would introduce bi-directional selectivity for the entire joint angle range (Fig. 2B). Note that the supplementary hairs of the extended bi-directional hair field were not linked to the position IN.

Neuron model for a high-pass filter

The mechanosensory neuron of each tactile hair has a phasic response at the moment when the hair reaches its maximum deflection, resulting in the spike rate peak shown in Figs. 3A,B. For fully deflected hairs ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi = {90}^{\circ }$$\end{document} ), the steady-state spike rate of the sensory neuron ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{\text {ss}}$$\end{document} ) is constant for all hairs at all times. In order to exploit the phasic behaviour of the sensory neuron, a single high-pass filter is integrated in series with the sensory neuron, whose cut-off frequency ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{c}$$\end{document} ) is just above \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{\text {ss}}$$\end{document} . Consequently, only the phasic response of the sensory neuron can trigger spikes from the high-pass filter. And thus, the neuron only fires when the joint angle exceeds the receptive field of the tactile hair and the hair reaches its maximum deflection. The joint velocity can then be extracted from the collective spike rate of the high-pass filters.

The LIF model, governed by Eqs. (9), (10), and (11), can operate as a high-pass filter (Mastella and Chicca 2021). Due to the leaky term in Eq. 10, the LIF only spikes if a specific input frequency is reached ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{c}$$\end{document} ). The cut-off frequency, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {c}$$\end{document} , depends on both \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} . By setting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} to a specific value and manually adjusting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {c}$$\end{document} can be set slightly above \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{\text {ss}}$$\end{document} . Filtered spikes from a hair field converge onto a single LIF neuron, which generates a spike in response to each input spike. This IN integrates the spike rates into a unified output, representing the velocity. A LIF generates a spike in response to a single input spike if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega> V_\text {T} - E_\text {L}$$\end{document} .

Results

Layer one: Hair field layer

Within the hair field layer, each joint was associated with two distinct hair fields, each containing a total of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h}$$\end{document} individual hairs. Within a hair field, each hair was sensitive to joint deflection in its own receptive field. Within this receptive field, the hair deflection angle was translated into an electrical current and transmitted to a single mechanosensory neuron designated to that hair. Subsequently, the mechanosensory neuron transformed an increment change in current into a phasic-tonic spiking response. This process mirrored the mechanosensitive dendrites found in actual tactile hair sensory neurons (Fig. 1B, Tuthill and Wilson 2016).

Replicating sensory spiking dynamics

The AdEx neuron was used to model the electrophysiological dynamics of tactile hair sensory neurons observed by Okada and Toh (2001), as depicted in Figs. 3A,B. This was achieved by replicating their experimental procedure and adjusting model parameters to match the observed spiking dynamics. Since the transient response of the AdEx model proved to decay faster than the experimental data, similar adaptation dynamics in the model required a relatively high \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _\omega$$\end{document} (approximately 600 ms). To prevent prolonged neuron suppression and maintain responsiveness to rapid changes in joint angle direction, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _\omega$$\end{document} in the AdEx neuron was reduced from 600 ms to 50 ms. A larger \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _\omega$$\end{document} caused the adaptation current to decay too slowly, suppressing phasic spiking upon reactivation following brief deactivation periods, such as those occurring during rapid reversals in movement direction. This would have mainly affected downstream velocity encoding, but not position encoding. Additionally, the corresponding angular velocities of hair deflection were quadrupled by dividing the total stimulus duration by four. This adjustment allowed us to replicate the general dynamics of the phasic-tonic response, albeit with different adaptation speeds. The modified procedures from the study of Okada and Toh (2001) were:

1. A ramp-and-hold function with linear increase of hair deflection at different speeds, while keeping the hold angle constant. Hairs were deflected from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${0}^{\circ }$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${37}^{\circ }$$\end{document} at five velocities. These were \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${980}\frac{\circ }{{s}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${604}\frac{\circ }{{s}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${352}\frac{\circ }{{s}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${188}\frac{\circ }{{s}}$$\end{document} , or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${96}\frac{\circ }{{s}}$$\end{document} .
2. A ramp-and-hold function linear increase of hair deflection at constant speed, but differing hold angles. Deflection velocity was \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${240.4}\frac{\circ }{{s}}$$\end{document} . The five hold angles used were: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${60}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${46}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${34}^{\circ }$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${23}^{\circ }$$\end{document} , or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${15}^{\circ }$$\end{document} . A grid search was employed to vary the values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}/{\phi }$$\end{document} and b within the ranges of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${10}\frac{{pA}}{\circ }$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${150}\frac{{pA}}{\circ }$$\end{document} and 50 pV to 350 pV respectively, with a total of 8 steps. The parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}/{\phi }$$\end{document} influences the current entering the neuron, thereby affecting the response strength. The parameter b modulates the degree of adaptation, impacting the relative strength of the phasic peak. Together, these parameters regulate both the phasic peak and steady-state frequencies, which are the key response metrics of the model. Several other parameters can be modified to achieve the desired spiking response (see supplementary Fig. S3, S4, and S5). The error metric mean absolute error (MAE) was computed to assess the deviation of the model response from the electrophysiological dynamics observed in the experimental study for critical spike rates (maximum and steady state), with:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {MAE} = \frac{1}{20}\sum _i^{10}(\mid \hat{x}_{\text {ss}}^{(i)} -x_{\text {ss}}^{(i)}\mid +\mid \hat{x}_{\text {max}}^{(i)} -x_{\text {max}}^{(i)}\mid ) \end{aligned}$$\end{document}

for the ith ramp function trial. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{\text {ss}}^{(i)}$$\end{document} refers to the spike rate at the steady state, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{\text {max}}^{(i)}$$\end{document} represents the maximum spike rate, and the notation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{x}^{(i)}$$\end{document} indicates values for the AdEx model while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x^{(i)}$$\end{document} indicates values extracted from Figs. 3A,B.Table 1. Model type and parameter values for the sensory neuron, position INs and velocity INs after optimization**IN**TypeC \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_\text {L}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_\text {L}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta _\text {T}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_\text {T}$$\end{document} a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _w$$\end{document} b \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} SensoryAdEx200 pF2 nS \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${-70}\hbox { mV}$$\end{document} 2 mV \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${-50}\hbox { mV}$$\end{document} 2 nS50 ms264 pV--PositionLIF-- \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${-70}\hbox { mV}$$\end{document} - \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${-50}\hbox { mV}$$\end{document} ---120 ms1 mVVelocityLIF-- \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${-70}\hbox { mV}$$\end{document} - \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${-50}\hbox { mV}$$\end{document} ---5 ms10.8 mVNote: This table presents the model parameters for different input types (sensory, position, velocity) after optimization

Optimization and performance

Before optimizing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}/{\phi }$$\end{document} and b, the initial AdEx parameters were adopted from Naud et al. (2008). In order to account for the high sensitivity of the mechanosensory neuron to small hair deflections, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_\text {L}$$\end{document} was reduced from 12 nS to 2 nS. This decreased the rheobase current, enabling the neuron to detect weak input currents. Additionally, the time constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _\omega$$\end{document} was reduced to 50 ms, as explained in the previous paragraph. Fig. 3E shows the MAE (computed using Eq. (12)) in each iteration of a parameter sweep over b and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}/{\phi }$$\end{document} . A minimum error of 20 was observed for several results, from which the values \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I/\phi = {50}\frac{{pA}}{\circ }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b = {264} \hbox { pV}$$\end{document} were randomly selected. The optimal parameters for the AdEx model are listed in Table 1.

Fig. 3. Phasic-tonic spike response of the real and modeled mechanosensory neuron. A. Spike response of mechanosensory afferents from the lateral scapal hair plate in the American cockroach (Okada and Toh 2001). The plot shows the spike rate over time for a tactile hair deflected by a ramp-and-hold function with a constant velocity of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${60.1} \frac{\circ }{{\textrm{s}}}$$\end{document} from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${0}^{\circ }$$\end{document} to various end angles. B. Response to ramp-and-hold deflection of a tactile hair with a constant hold angle of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${37}^{\circ }$$\end{document} , but varying angular velocities. C. Stimulus protocol (top) and corresponding spike responses of the AdEx model (bottom). The stimulus protocol is the same as for A, except with four-fold ramp velocity and shorter hold time. D. Stimulus protocol (top) and corresponding spike responses of the AdEx model (bottom). The stimulus protocol is the same as for B, but with four-fold ramp velocities and shorter hold time. E. Heatmap illustrating the mean absolute error (MAE) between significant spike rates (peak and steady-state) of the modeled and experimental sensory neuron for a parameter sweep across \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}/{\phi }$$\end{document} and b. The MAE was computed according to Eq. (12)

The results of both hair deflection experiments are illustrated in Figs. 3A,C (varying the hold angle) and 3B,D (varying angular velocity of the ramp) using the optimal parameter set. The dynamics can be directly compared to the experimental data on hair field afferents of the lateral scapal hair plate of the American cockroach (Figs. 3A,B, Okada and Toh 2001). The stimulus protocols employed for the experimental and model data were comparable, although the former incorporated four-fold angular velocities and a shorter hold time. Despite these adjustments, the modelled sensory neuron exhibited phasic-tonic dynamics similar to that observed in the experimental data. During the ramp phase of the stimulus, the spike rate showed a linear increase, reaching a peak as the ramp reached the constant hold angle. Subsequent to the peak, the spike rate converged to a steady state. As the stimulus was ramped back to zero, the spike rate showed a linear decrease until it reached zero. In contrast to the considerable noise present in the experimental data, the model is deterministic and, therefore, devoid of noise. Furthermore, the model reached a steady state within 100 ms, whereas the experimental data displayed a continuous decay throughout the entire hold phase.

In trials with varying hold angles (Fig. 3A,C) the transient peaks reached approximately 110% of the steady-state conditions in the model. In contrast, the peaks in the experimental data were approximately 150% of the steady-state spike rate. Moreover, in the experimental data, the low-end deflection angles (e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta = {15}^{\circ }$$\end{document} ) only produced a response during the onset of the ramp function, not during the steady state. This behavior was not observed in the simulated model, where low deflection angles induced a steady-state spike response. Except for very small hold angles, the steady-state spike frequencies of the model aligned closely with experimental results.

Trials with varying ramp velocities (Figs. 3B,D) exhibited a higher degree of similarity. The steady-state spike rate was slightly higher than the noise range of the experimental data. Additionally, the peaks were in a slightly wider range: 108 Hz - 199 Hz for the model compared to 118 Hz - 166 Hz for the experiment.

Layer two: Position interneurons (INs)

This section commences with an evaluation of the characteristics of the bi-directional hair field, which is followed by an analysis of the performance of the position INs. Two position INs were integrated in each joint, with the function of integrating all spikes originating from a given hair field into a position-dependent spike train. These spike trains were able to encode the varying joint angles over time through their spike rate.

Bi-directional hair field

In each hair plate, 100 hairs were arranged to form two hair fields, each comprising 50 hairs ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{\text {h}} = 50$$\end{document} ). For each joint, the bounding joint angles of the bi-directional hair field, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _\text {j}^{(max)}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^{\text {(min)}}_{j}$$\end{document} , were set as the maximum and minimum values of the joint’s working-range. The overlap was defined as follows: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta ^{\text {(olhf)}} = \theta ^\text {(ol)} = {0.1}^{\circ }$$\end{document} . As was clarified in the methods section, half of the hairs in a hair field were connected to position INs. Accordingly, Fig. 4A includes only hairs with synaptic connections to position INs. The time-varying plot displays a 5-second time course for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} joint angle of the right front leg, R1. The dotted line represents the resting angle, situated in the middle of the joint range. Red and blue dot sequences form a hair field raster plot, where each dot denotes a spike event for the ith hair, respectively. Due to the design of the hair field, proprioceptive hairs numbered 1 to 25 are sensitive to joint angles in the negative domain, while hairs 51 to 75 are sensitive to joint angles in the positive domain. Due to the large number of spikes, dots representing individual spike events coalesce into a continuous line.

Fig. 4. Position IN results. A. Raster plot displaying the response of a hair plate ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h} = 50$$\end{document} ) to a 5-second movement sequence of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} joint (thorax-coxa) of the right front leg (R1). Spike events from proprioceptive hairs in the anterior (blue dots) and posterior (red dots) parts of the joint’s working range encode increasingly protracted or retracted postures, respectively. The high spike density results in the dots merging into a continuous line. B. Heat map depicting the mean squared error (MSE) between experimentally obtained joint angles and model predictions across a parameter sweep of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} . The results are averaged over 78 trials and 18 joints, totaling 1404 joint angle time courses. C. Spike rate time courses of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pos_+$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pos_-$$\end{document} INs in response to the same joint angle movement as shown in panel A. D. Combined and z-normalized spike rates of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pos_+$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pos_-$$\end{document} INs from panels A and C closely follow the joint angle time course. The model response and experimental data are plotted in arbitrary units due to z-normalization.

Fig. 4A highlights the bi-directional sensitivity of the hair field, with mirrored encoding observed in the positive and negative domains. In the positive domain, increasing joint angles prompted additional neurons to activate, whereas in the negative domain, the opposite was true. At the resting angle, only the initial sensory neurons (25 and 26) activated due to the pre-defined overlap. Ultimately, the close alignment between joint angle and hair field activity was evident.

Optimization and performance

To assess the performance of the position INs, we interpolated the experimental joint angles to match the timestep of the model ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$dt = {5} \hbox { ms} \rightarrow {0.25} \hbox { ms}$$\end{document} ). Subsequently, we subtracted the spike rate of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pos_-$$\end{document} IN from the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pos_+$$\end{document} IN and then z-normalized the resulting time series ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\psi }_{jk}(t_\text {i})$$\end{document} ). For direct comparison, we z-normalized the corresponding joint angle time course ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\psi _{jk}(t_\text {i})$$\end{document} ), as illustrated in Fig. 4D. The squared difference between these z-normalized time courses was then averaged over each discrete time step i, joint angle j and trial k, yielding the MSE for one joint angle time course:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {MSE} = \sum ^{78}_{k=1} \sum ^{18}_{j=1} \sum ^{N_\text {steps}}_{i=1} \frac{(\hat{\psi }_{jk}(t_\text {i}) - \psi _{jk}(t_\text {i}))^2}{18\times 78\times N_\text {steps}} \end{aligned}$$\end{document}

The MSE, averaged over 78 trials and 12 joint angles, was computed at each iteration of a parameter sweep for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} , and the results were visualized in a heat map (Fig. 4B). The optimal parameters, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau = {120} \hbox { ms}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = {1} \hbox { mV}$$\end{document} yielded an MSE of 0.03053. A decrease in MSE was observed while lowering \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} . The synaptic weight was not reduced further because the maximum observed spike rate (Fig. 4C, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{max}$$\end{document} = 322 Hz) was already low. Further reduction would have introduced more stair stepping. The time constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} influenced the MSE little at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} values of 2–10 mV. While for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = {1} \hbox { m V}$$\end{document} , low \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} values increased the MSE. At these values, the spike rate was too low to allow for detailed encoding by spike rate modulation, due to pronounced stair stepping. The complete optimal parameter set of the position IN is listed in Table 1. Fig. 4C illustrates a 3-second joint angle time course and spike rates for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} joint of R1. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pos_+$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pos_-$$\end{document} neurons fired proportionally to the joint angle relative to the resting angle. Spike rates fluctuated between 0 and 322 Hz and followed the joint angle time courses (though inverted in case of the hairs in the posterior working-range). Note that the degree to which fine details of the joint angle time course were lost in spike rate modulation depended, apart from the model parameters mentioned above, on the number of proprioceptive sensory neurons per hair field. If spike rates in the kHz range are required for further processing, as is the case for the subsequent layer in the companion paper (companion paper, van der Veen et al. 2025), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} can be increased (e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} = 25mV for the companion paper). To visually assess joint angle encoding, the response of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pos_-$$\end{document} was subtracted from the response of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pos_+$$\end{document} and then z-normalized. This time course is plotted alongside a z-normalized joint angle time course in Fig. 4D. Visually, the spike rate closely resembled the joint angle time course. However, the model response ‘overshot’ the joint angle time course at moments of sudden changes in joint angle (e.g., t = 0.5 s, 1.5 s, 2.5 s). Additionally, there was some stair stepping in the spike rate.

To test whether position IN performance was the same irrespective of leg or joint types, a two-way Analysis of Variance (ANOVA) (Fisher 1970) was conducted to examine the effect of the factors ‘joint type’ ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} ) and ‘leg type’ (front, middle, and hind) on the MSE. The results showed no significant main effect of factor ‘leg type’ (F(2, 69) = 2.52, p = 0.084). However, there was a significant main effect of factor ‘joint type’ (F(2, 69) = 111.62, p = 9.62 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-32}$$\end{document} ), indicating that position encoding varied significantly among joints. A post hoc comparison using the Tukey’s honestly significant difference (HSD) test (Tukey 1949) indicated that the MSE of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} joint was significantly different from the MSE of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} joints. Table 2 lists the specific MSE values for the three leg types and three joint types. The table reflects that position encoding was statistically indistinguishable between leg types, but significantly worse for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} joint compared to the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} joints. Additionally, a significant interaction between leg types and joint types was observed (F(2, 69) = 11.11, p = 4.90 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-8}$$\end{document} ). A post hoc comparison using the Tukey’s HSD test indicated that for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} joint, the front legs performed significantly worse, and for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} joint, the middle legs performed significantly worse.Table 2MSE and standard deviation (SD) for the position INs, grouped by factors legs and joints of the two-way ANOVAJointMSESDfront0.03170.0071middle0.03100.0105hind0.02960.0060 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} 0.02660.0046 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} 0.03890.0075 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document} 0.02680.0047

Layer two: Velocity interneurons (INs)

The objective of this section is to test how a velocity-sensitive first-order IN could encode joint angle velocity independent of joint position, despite its presynaptic input elements being position-sensitive. As with position INs, the instantaneous joint angle velocity could be encoded through spike rate modulation. With regard to the encoding properties of the mechanosensory neurons, it could be observed that movement was signaled by transition of joint angle from one hair’s receptive field to the next. In order to encode such transition events in a systematic manner, we decided that the phasic response of the sensory neurons should be leveraged. The phasic response component of a hair field afferent was driven by the initial deflection of the hair, before the adaptive current of the AdEx model built up. It was possible to integrate several phasic response components within an entire hair field in a first-order IN by emphasizing subsequent joint angle transitions from one receptive field to the next. As this mechanism strongly depended on the phasic response property of the AdEx, its time constant needed to be sufficiently short which, in turn, required the temporal fourfold scaling of the afferent response (see Fig. 3 and associated text).

Optimization and performance

First, parameters were optimized to set the high-pass filter properties of the velocity IN such that the cut-off frequency \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {c}$$\end{document} was just above the steady state frequency \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {ss}$$\end{document} . Since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {c}$$\end{document} depends on both \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} was arbitrarily set to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${5} \hbox { ms}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} was increased manually until \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {c}$$\end{document} was reduced to just above \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {ss}$$\end{document} , resulting in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} = 10.8mV. With this value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} , a minimal phasic peak in the AdEx neuron would be sufficient to excite a spike in the high-pass filter, making the velocity IN sensitive to low velocity.

Fig. 5A illustrates a linear relationship between the spike rate of a velocity IN and joint angular velocity. It is evident for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} = 10.8mV, small angles are encoded while for low \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} (9mV), angular velocity was not or weakly encoded at 100 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\circ }{\textrm{s}}$$\end{document} and below. Increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h}$$\end{document} (from 25 to 50) resulted in a higher overall spike rate, stronger linearity (higher \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^2$$\end{document} value). Considering most accurate, linear encoding of velocity, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h} = {50}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} = 10.8mV were selected as the optimal parameter set. All model parameters of the optimal velocity INs are listed in Table 1. For these parameters, the spike rate increased with a slope of 2.02 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{\circ }$$\end{document} . The linear velocity dependency was very much like that velocity-sensitive descending INs iONv and cONv of the stick insect antennal mechanosensory pathway (Fig. 4 in Ache et al. 2015), except for the lacking resting activity.

Fig. 5. Velocity IN results. A. For velocity INs, the output spike rate increases linearly with angular velocity. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h}$$\end{document} were varied to test their impact on the spike rate. The four variants shown here differ in threshold velocity for reliable encoding and the slope of the linear relationship. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r^2$$\end{document} value of a linear regression fit is given for each variant. B. Spike events superimposed on the joint angle time course reveal consistent encoding of positive slope by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_+$$\end{document} INs (blue) and of negative slope by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_-$$\end{document} INs (red). C. The spike rate of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_-$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_+$$\end{document} INs (red and blue, respectively) in response to a concurrent change in joint angle velocity (black). D. Combined and z-normalized spike rate of both the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_+$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_-$$\end{document} INs in response to a representative experimental angular velocity record. Both the model response and experimental data are plotted in arbitrary units due to z-normalization. A persistent time lag of approximately 0.025 s is observed in the model response.

Fig. 5B plots individual spikes for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_-$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_+$$\end{document} INs, overlaid on the joint angle time course for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} joint of the right front leg, R1. Evidently, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_-$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_+$$\end{document} neurons spike during opposite movement directions. The figure illustrates a delayed offset after sudden changes in joint angle (peak or trough). During these changes, the current in the mechanosensory neuron can decrease but remain high. If adaptation has not fully occurred, the sensory neuron spike rate can still exceed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {c}$$\end{document} . Consequently, a velocity IN may continue to fire for a short duration even after the joint angle has changed direction, leading to the observed delayed offset and reduced accuracy. Additionally, since velocity IN spikes only occur near the edges of the receptive field, there is no sensitivity to movements within the receptive field, further lowering accuracy.

Fig. 5C plots the angular velocity over time in conjunction with the spike rate for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_+$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_-$$\end{document} INs. Fig. 5D plots the combined and z-normalized spike rates of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_+$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_-$$\end{document} INs. The model response lags behind the experimental data by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t \approx {0.025} \hbox { s}$$\end{document} . Since the error measures of the position and velocity INs are similar, their values can be directly compared. Overall, velocity INs (MSE \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$= {0.3421}$$\end{document} ) were less accurate than position INs (MSE \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$= {0.03053}$$\end{document} ). If shifted by 0.025 s, to correct for the lag, the MSE improved to 0.0910, which was still by 0.0605 higher than in position INs.

To further assess the accuracy of the model, the velocity INs were treated as binary classifiers. Spikes of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_+$$\end{document} type IN during positive and negative movements (i.e., increasing and decreasing joint angles) were denoted as true positive (TP) and false positive (FP), respectively. Similarly, spikes of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_-$$\end{document} type IN during positive and negative movements were labeled as false negative (FN) and true negative (TN), respectively. Since the number of positive (P) and negative (N) occurrences was expected to be balanced, a simple accuracy calculation was sufficient:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} TPR = \frac{TP}{TP+FN} \\ TNR = \frac{TN}{TN+FP} \\ ACC = \frac{TPR + TNR}{2} \end{aligned} \end{aligned}$$\end{document}

where the true positive rate (TPR) is the sensitivity and the true negative rate (TNR) is the specificity. An accuracy of zero indicates that all spikes occur at incorrect times, while an accuracy of one implies that all spikes occur at the correct time. The averaged results over 18 joint angles and 78 trials are provided in Table 3. The accuracy is high (0.914), and both TPR and TNR are balanced, indicating that both neurons ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_+$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$vel_-$$\end{document} ) perform similarly accurate.Table 3. Confusion matrix and accuracy metrics for the velocity INs if regarded as a binary classifierPredictedActualPositiveNegativePositive1627811155707Negative1551671689086AccuracyTPRTNR0.9140.9130.916

Discussion

Our results demonstrate the use of computationally compact integrate-and-fire neuron models for distributed encoding of joint kinematics. Our focus on proprioceptive hair fields as mechanosensory inputs does not limit the generality of the model for distributed proprioception in insects because two major computational features are common to most insect proprioceptors: (i) the spike rate adaptation in phasic-tonic spike responses and the associated fractional encoding of the sensory magnitude and its derivative; and (ii) range fractionation of the input magnitude in an array of strictly local sensory neurons. The first layer of our model captures both of these features (Figs. 2, 3). The second layer established linear encoding of position (Fig. 4) and velocity (Fig. 5) in first-order INs, despite mixed position-velocity encoding of their presynaptic inputs.

Layer one: Hair field

Proprioceptive hair fields come in two types, either as two-dimensional patches (hair plates) or linear rows of hairs (hair rows), sometimes with both types occurring on the same limb segment (e.g. stick insect antennae: Krause et al. 2013; stick insect legs: Wendler 1964; Fig. 1A). As hair fields are often arranged in pairs at opposite sides of the same limb segment, the sensory array of our hair field layer was designed as a pair of opposing hair rows. Accordingly, the bi-directional hair field model in Fig. 2 is a structural template of real hair fields in that the angular working-range of the joint is sampled by two complementary, idealised sensory arrays with equal spacing of hairs.

Structural idealisations of our sensory array concern, for example, the uniform length of tactile hairs, neglecting known length variation (Pringle 1938) and functional subdivisions of two-part hair fields, as in the trochanteral hair field of the stick insect (Schmitz 1985b, 1986). Secondly, a linear relation between hair angle and joint angle was assumed, leading to a linear dependence of tonic spike rate and joint angle, and neglecting non-linearities such as the threshold deflection angle for tonic spike activity (see Fig. 3A and below). Finally, our assumption of a linear sensory array (hair row) neglects overlapping receptive fields for multiple hairs that occur for two-dimensional hair plates. Computationally, it is easy to expand our structural template so as to include non-equal spacing and/or multi-hair overlap. On the other hand, there is evidence for tactile hairs being arranged at regular intervals (Tuthill and Wilson 2016; Hannah-Alava 1958).

The phasic-tonic dynamics of a position-dependent response with spike rate adaptation was modelled after experimental data from afferents of the lateral scapal hair plate of the American cockroach (Figs. 3A,B; Okada and Toh 2001). Other than earlier studies that modelled afferent hair field spikes with a uniform random number generator (Ache and Dürr 2015) or Poisson process (Gollin and Dürr 2018), we employed an integrate-and-fire neuron model. This model family captures the non-linear transformation that occurs between the stimulus (e.g., deflection of the hair) and the neuron’s response, including sub-threshold temporal integration and a spike threshold that is set by a biologically meaningful parameter rather than by mathematical parameters (as in a Poisson process). The phasic-tonic properties associated with spike rate adaptation were incorporated through the AdEx model (Figs. 3C, D). In contrast to non-spiking models of non-linear afferent dynamics (Szczecinski et al. 2021), the AdEx model directly generates spikes rather than capturing only changes in spike rate. Moreover, it accepts raw input of joint angles without undergoing pre-processing via high-pass and low-pass filters, as in (Ache and Dürr 2015). Here, the sole pre-processing step involved is converting joint angles into hair deflection angles, subsequently transforming them into an input current.

The adaptation term in the AdEx model allowed for a combined phasic-tonic response and therefore copied the essential features of the observed experimental spike dynamics, though in a mathematically compact manner. As a trade-off compared to more physiological Hodgkin-Huxley-type models of mechanotransduction (e.g., for spider slit sensilla; French et al. 2002) integrate-and-fire models are not suitable for sensitivity analyses addressing the function of particular ion channels.

The transient time course of spike rate adaptation in the physiological recordings from cockroaches (Figs. 3A,B) proved to be slower than what could be achieved by the AdEx model. The best AdEx model fit required a long time constant ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau \approx {600} \hbox { ms}$$\end{document} ). Apart from leading to significant delays in subsequent layers, it mainly affected the tuning of the downstream velocity INs. After a sudden change in joint angle, the mechanosensory neuron remained suppressed due to adaptation for a relatively long time which, in turn, lead to inconsistent performance in the velocity INs, and reduced accuracy in position INs. To alleviate this effect, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document} was adjusted to 50 ms, enhancing adaptation speed while preserving phasic-tonic behaviour. After optimization, there was a non-negligible MAE of 20 Hz. Moreover, our hair field afferents (Fig. 3A) encoded considerably lower deflection angles than observed for the cockroach hair field (Fig. 3C). Finally, our model is deterministic in that it does not include any noise.

By choosing to optimize the spike dynamics according to a proprioceptive hair plate of the American cockroach, we assumed that these dynamics would be consistent across different species and limbs. However, given the fact that even hairs within the same hair plate may differ with regard to details of their spike response (Wong and Pearson 1976) a general-purpose model needs to be tunable. Indeed, the AdEx model allows simple tuning of phasic peak and tonic steady-state spike rates by adjusting specific model parameters. To demonstrate this, the supplementary material shows how to adjust phasic and/or tonic responses as desired. For instance, varying a modifies \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {ss}$$\end{document} while only minimally affecting the phasic peak (Fig. S3). Increasing C while decreasing a alters the strength of the phasic peak without changing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {ss}$$\end{document} (Fig. S4). Additionally, the AdEx model is capable of generating a purely tonic response by setting b to zero and modifying C (Fig. S5). Since most proprioceptors display a phasic and/or tonic response, this demonstrates the model’s versatility in simulating a wide range of proprioceptors and their dynamics.

We conclude that the AdEx model may serve a general template for spiking proprioceptor afferents with phasic-tonic response dynamics. To obtain an even more biologically accurate model, future studies could explore methods to account for variations in hair length, analyze possible nonlinear relationships between joint and hair angle, or optimize hair distribution. The latter could reduce redundancy and/or entail the implementation of the efficient coding hypothesis (Barlow et al. 1961). This hypothesis states that sensory systems are configured to efficiently encode the dynamic sensory stimuli that the organism encounters in its environment. Therefore, these systems are unlikely to adhere to a strict linear structure; instead, they are likely arranged to maximize coding efficiency for the most common sets of stimuli encountered by the organism. This concept can also be applied to the organization of the hair fields, suggesting that hairs can exhibit variation in receptive field size based on the probability distribution of joint positions.

Layer two: Position interneurons (INs)

Position-sensitive second-order INs were modelled to extract tonic proprioceptive information from the hair field layer. The design drew inspiration from so-called ‘simple position-sensitive’ descending INs of the stick insect antennal system (Ache and Dürr 2013, 2015). Together with ‘dynamic position-sensitive’ descending INs that spike whenever the antenna is moving across a particular joint angle range, they are well-documented examples of position-encoding first-order INs in general. Position-dependent encoding of joint angles of foot position has also been found in thoracic INs of the stick insect (Brunn and Dean 1994; Berg 2014) but were mostly characterised with cyclic leg movement stimuli that combined position- and motion cues. Concerning the putative proprioceptors driving these IN responses, the trochanteral hair field has been shown to be the sole sensory input to the feedback control loop at the coxa-trochanter joint (Schmitz 1986). Accordingly, ablations of this hair plate affects leg movement in different species, including swing height of walking Carausius morosus (Theunissen and Dürr 2013), and the amplitude of searching movements in Medauroidea extradentata (Berg et al. 2013).

Other than an earlier model by Ache and Dürr (2015), the present study uses a spiking neuron model to implement position-encoding first-order INs. Consistent with our study, they found that employing two hair fields connected to two separate INs allowed to estimate joint angle reliably from a natural stimulus.

The non-zero MSE in the position layer ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$MSE = {0.03053}$$\end{document} ) arose due to two main factors: Overshoot (Fig. 4D, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t = {0.5},~{1.5},~{2.5}$$\end{document} ) and stair-stepping (Figs. 4C,D). Overshoot occurred due to the phasic peak in the sensory neuron activity during changes in joint angle direction. The position INs adopted this peak and, therefore, overshot relative to the experimental joint angle. The overshoot was exacerbated by increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} , as evidenced by the rise in MSE depicted in Fig. 3B. To mitigate this overshoot, one could introduce a low-pass filter between each sensory neuron and its corresponding position IN. Stair-stepping arose as an artefact from the simulation time step and the calculation of the spike rates. This artefact can be reduced by shortening the time step or by smoothing the spike rate time course.

Layer two: Velocity interneurons (INs)

As a complement to position INs, velocity INs were designed to extract the phasic component of afferent activity, and to spike only if the joint angle changed. Velocity-encoding spiking interneurons in stick insects were first recorded by Büschges (1989). The highly linear velocity INs from Fig. 5 were inspired by ‘ON-type velocity sensitive’ INs in the stick insect antenna (Ache and Dürr 2013; Ache et al. 2015). At least one of these receives proprioceptive input from pedicellar hair fields (Jaske et al. 2021). There is also evidence for velocity-sensitive INs in the antennal systems of crickets (Gebhardt and Honegger 2001) and the fruitfly (Mamiya et al. 2018). In comparison with an earlier model (Ache and Dürr 2015), our present model replaces filter blocks by a SNN. As their physiological counterparts, our velocity INs demonstrated linearity between spike rate and angular velocity (Fig. 5A) and precise encoding of angular velocity in response to realistic stimuli (Fig. 5D, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$MSE = {0.3421}$$\end{document} ). Velocity encoding was direction-selective (Figs. 5B,C) with a high accuracy of 0.914.

Strong linearity was achieved only when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} equalled \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${10.8} \hbox { mV}$$\end{document} . At this synaptic strength, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {c}$$\end{document} marginally exceeded \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {ss}$$\end{document} , allowing spikes from small phasic fluctuations to transmit through the high-pass filter. However, for lower values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {c}$$\end{document} increased, resulting in complete blockage of spikes from small phasic fluctuations by the high-pass filter and zero postsynaptic spikes at low angular velocities (Fig. 5A). Moreover, increasing the number of hairs per hair field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h}$$\end{document} also enhanced encoding of low angular velocities. Doubling \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h}$$\end{document} (from 25 to 50) halved the receptive field size and doubled the rate of change of hair angles for a given change in joint angle. This amplified the phasic peak for low angular velocities. Thus, the optimal parameter set would be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = {10.8} \hbox { mV}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h} = \infty$$\end{document} . However, due to a biological constraint on the number of hairs, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h} = {50}$$\end{document} was considered the maximum number of hairs per hair field (Pringle 1938).

The apex of the phasic peak occurred when a hair reached maximum deflection (Fig. 3D). Consequently, the second half of the phasic peak lagged maximum hair deflection. This resulted in a slight delay ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx {0.025} \hbox { s}$$\end{document} ) in the velocity INs’ spike rate compared to the actual joint velocity. This delay was clearly evident in Fig. 5B and D, and largely contributed to the accuracy of 0.914. Upon correction of this delay, by manually shifting the model response by 0.025 s, the MSE decreased from 0.2812 to 0.0910 (compared to 0.03053 for the position INs). This implies that the velocity INs were less accurate than the position INs. Position INs were capable of detecting joint angle changes within the receptive field of a hair, whereas velocity INs only responded to changes from one receptive field to another. Consequently, the sensitivity of velocity INs was directly linked to the number of receptive field edges and therefore \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_\text {h}$$\end{document} .

There is some discrepancy between the ON-type neurons found in the stick insect and our velocity IN implementation. Firstly, both hair fields had to be extended (Fig. 2) to achieve sensitivity across the entire joint angle working range. Secondly, our velocity INs were directionally selective, unlike the ON-type neuron (Ache et al. 2015). Ache and Dürr ( 2015) had already suggested two ways by which a velocity IN could lose the directional selectivity of its afferent input (post-excitatory rebound; additional interneurons). A further option would be to fully connect all hairs from a extended bi-directional hair field (Fig. 2B). Here, we renounced any of these options because direction-selectivity demonstrated distinct advantages in the subsequent network layers (see companion paper: van der Veen et al. 2025). Thirdly, the current network architecture solely relies on afferent input from hair fields. However, position and velocity information at the insct femur-tibia joint is supplied by phasic-tonic afferent responses of the femoral chordotonal organ (e.g., Zill 1985; Hofmann et al. 1985; reviewed by Field and Matheson 1998). It could be interesting to explore whether incorporating additional proprioceptors would mitigate the issues raised.

Future work

We conclude that two fundamental aspects of distributed proprioceptive encoding of joint kinematics can be modelled by a two-layered SNN: phasic-tonic encoding and the extraction of joint angles and joint angle velocity at high accuracy. Future research may expand on these results, e.g., by tuning the AdEx neuron model parameters to different proprioceptor types (e.g., scolopidial neurons of chordotonal organs; Tuthill and Wilson 2016; Mamiya et al. 2018), by introducing multimodal integration among different proprioceptors (Gebehart et al. 2022), or by detailed models of distinct classes of first-order interneurons (Ache and Dürr 2013, 2015; Berg 2014).

The companion paper (van der Veen et al. 2025) expands the present network by incorporating ‘movement primitive’ neurons that are designed to extract more complex intra-leg information about step cycle phases and/or transitions between them. Additionally, it introduces posture neurons that are capable of encoding higher-order information about whole-body posture from movement primitive neurons of all six legs. In combination with the present study, the proposed four-layered SNN links aspects of primary proprioception, inter-segmental proprioceptive pathways, and internal representation of whole-body posture.

Supplementary Information

Below is the link to the electronic supplementary material.

Supplementary Material 1: The AdEx neuron dynamics.The AdEx model response to a sustained current, I, governed by Eqs. (6, 7, 8). This current is integrated into the membrane voltage, V. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V = V_\text {T} = {-50} \hbox { mV}$$\end{document} , the membrane voltage spikes and resets to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_\text {L} = {-70} \hbox { mV}$$\end{document} . During a spike event, 0.264 nA is added to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} , counteracting the input current. The model exhibits lower spike frequencies with increasing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega$$\end{document} until reaching an equilibrium spike rate. The model parameters are given in Table 1 (sensory neuron).

Supplementary Material 2: The LIF model dynamics. The LIF model responds to a constant spike rate (333.3 Hz), governed by Eqs. (9, 10, 11). A presynaptic spike, Pre, increases the membrane potential by 10 mV, and if no spikes are present, the membrane voltage, V, decays back to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_\text {L} = {-70} \hbox { mV}$$\end{document} . When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V> V_\text {T} = {-50} \hbox { mV}$$\end{document} , a postsynaptic spike, Post, is recorded, and V resets to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_\text {L} = {-70} \hbox { mV}$$\end{document} . The postsynaptic spike rate remains constant in response to a constant presynaptic spike rate, since the LIF model has no adaptation. At the shown input spike rate, every third presynaptic spike triggers a postsynaptic spike. The model parameters are given in Table 1 (velocity neuron)

Supplementary Material 3: The AdEx model response to varying a. The AdEx model’s spike response to a ramp-and-hold hair deflection at an angular velocity of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\circ }{\textrm{s}}$$\end{document} . The parameters for the AdEx model were taken from Table 1. The parameter a was varied to demonstrate that the tonic steady-state frequency can be reduced while maintaining a strong phasic peak. A purely phasic response can be achieved by combining the AdEx model with a high-pass filter in series.

Supplementary Material 4: The AdEx model response to varying C and a. The AdEx model’s spike response to a ramp-and-hold hair deflection at an angular velocity of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\circ }{\textrm{s}}$$\end{document} . The parameters for the AdEx model were taken from Table 1. The parameters C and a were varied to demonstrate that the strength of the phasic peak can be adjusted without altering \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {ss}$$\end{document}

Supplementary Material 5: The AdEx model response to varying C. The AdEx model’s spike response to a ramp-and-hold hair deflection at an angular velocity of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\circ }{\textrm{s}}$$\end{document} . The parameters for the AdEx model were taken from Table 1, with b set to zero to eliminate adaptation and yield a purely tonic response. The capacitance C was adjusted to achieve different values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\text {ss}$$\end{document} .

Supplementary Material 6

Supplementary Material 7