A Neuro-Musculo-Skeletal Model for Insects With Data-driven Optimization

Simulating the locomotion of insects is beneficial to many areas such as experimental biology, computer animation and robotics. This work proposes a neuro-musculo-skeletal model, which integrates the biological inspirations from real insects and reproduces the gait pattern on virtual insects. The neural system is a network of spiking neurons, whose spiking patterns are controlled by the input currents. The spiking pattern provides a uniform representation of sensory information, high-level commands and control strategy. The muscle models are designed following the characteristic Hill-type muscle with customized force-length and force-velocity relationships. The model parameters, including both the neural and muscular components, are optimized via an approach of evolutionary optimization, with the data captured from real insects. The results show that the simulated gait pattern, including joint trajectories, matches the experimental data collected from real ants walking in the free mode. The simulated character is capable of moving at different directions and traversing uneven terrains.

in real insects and in charge of their gait pattern. Researchers also modelled and animated virtual myriapoda by developing a hybrid animation system that combines kinematic and dynamic simulation 13 . While the body was driven by both rigid body simulation (for hard skeleton) and finite element method (for soft interlinks), the leg movements were synthesized in a kinematic fashion.
Robotics researchers used the simulation framework to speed up the process of robot design and control. A representative hexapod robot, RHex 14 , was designed to perform the double-tripod gait similar to the gait pattern observed in hexapod insects. This gait pattern allowed inherent locomotion stability. The open-loop controller was further improved by the introduction of Central Pattern Generator (CPG). The biological evidences show that the CPG is responsible for producing the basic gait pattern of real animals 15 . Bio-inspired robots, such as walking salamander 16 , swimming fish 17 and flying bird 18 , are controlled by mimicking the machinery of CPG. More details can be found in the review paper on applications of CPG for locomotion control in animals and robots 19 . However, existing models of CPG are generally constructed as a network of nonlinear oscillators, rather than spiking neurons. The advantage of using nonlinear oscillator is the inherent convergence to limit cycle or point attraction for most oscillator models, thus capable of endogenously generating the gait pattern without the high-level commands. However, the challenge from oscillator-based networks is how to connect the sensory information/high-level commands with the low-level control strategy. Previous method 11 used a precomputed look-up table. This inevitably creates a large discretized look-up table. In comparison, real organisms, including insects, use spiking neurons for the general purposes of environmental sensory, information transmission, decision-making and muscle actuation. Such observation inspires this work to explore the solution of using spiking neuron, as the alternative building block to nonlinear oscillator, to construct the neural system. Building a simulation model to match its real equivalence is a non-trivial task. On one hand, a complex model may introduce too many parameters, which could be impossible to tune by hand. On the other hand, a simplified model may not capture the important details of the model and thus reduce the credibility and capability of the virtual model. In this work, we propose a neuro-musculo-skeletal model for virtual insects (using specifically the ant as an example) and an automatic solution to the problem of parameter optimization. More specifically: • We propose a uniform encoding of sensory information, high-level commands and control strategy, with the spiking patterns of the neural system. Such an encoding paradigm provides the equivalent representation as real organism and shows the potential of using spiking neurons as the building block for accomplishing general tasks. The network also provides a continuous mapping from the sensory information/high-level commands to low-level control, in contrast to the discretized mapping in previous work 11 . • We use the evolutionary algorithm to automatically find the optimal values for both the controller and muscle parameters. Given the high dimensions of the model state spaces, we resolve the challenge by dividing the optimization into a two-stage procedure. The first stage optimizes the parameters of controller and muscles for the standard gait of walking straight forward, while the second stage progressively extends the skill repertoire by considering curve trajectories and uneven terrain. • We integrate the prior information (foot trajectories) of real ants during the optimization process via the evolution algorithm. The introduction of ground-truth data from real ants allows the quantitative comparison between the proposed model in this work and existing models 11 . The result shows that the trajectories generated from this model presents smaller deviations from the ground-truth data.

Method Overview
The simulation framework is composed of three parts: the spiking-neuron controller, the muscle actuation and the skeleton. The spiking-neuron controller receives the external inputs (sensory information and high-level commands) and sends out a series of spiking trains during each motion cycle and activates the muscle contraction. Each skeleton joint is connected with a pair of antagonistic muscles, generating the appropriate torques and resulting in the forward movements of the virtual insect. The parameters of the neuron controller and muscle actuation are updated with the technique of evolutionary algorithm, with the captured data of real insects. The overall algorithm of our method is presented in Fig. 1(a). Data availability. All data generated or analysed during this study are included in this published article (and its Supplementary Information files).

Neuro-Musculo-Skeletal Model
Neuro-controller Model. In real insects, spiking pattern of the neural system are used for the general purposes of information transmission and processing. The controller in the proposed model is modelled as a network of spiking neurons. The unit model is based on the model proposed in the work by Izhikevich 20 . The model is endowed with both biologically plausibility of Hodgkin-Huxley-type dynamics and computational efficiency of integrate-and-fire neurons. The equation for the spiking neuron model is: with the auxiliary after-spike resetting: v represents the membrane potential of the neuron. u represents the activation of K + ionic currents and inactivation of Na + ionic currents, providing negative feedback to v. This model has four dimensionless parameters: a, b, c, d. The model is capable of generating multiple spiking patterns, including regular spiking, intrinsically spiking, chattering and so on. We here choose the most common type of spiking pattern (regular spiking) and its corresponding parameters are: a = 0.02, b = 0.2, c =−65 mV, d = 8. One novelty of this work is to provide a uniform encoding of both sensory information, high-level commands and control strategy. We choose specifically the slope angle as the example of the sensory information and the action of turning as the example of high-level commands. For the slope angle θ s or the target orientation for next stride θ t , the input currents I are defined as The input currents I sensor , I target are fed into an input layer of spiking neurons. The input layer is connected with a central layer of 16 spiking neurons, which are then fully-connected to the motor neuron. The output from the central layer is added with a bias component (denoted in blue in Fig. 1(b)), which generates the standard gait pattern of walking straight forward. This constructs a neural network with spiking neurons ( Fig. 1(b)). The synaptic connection weights between the neuron layers determine how the external information adjusts the control signals from the bias neurons and eventually modifies the standard gait pattern.
For the control strategy, each muscle is activated by a unit neuron model, which produces the overall spiking pattern of the neuron group activating the single muscle. Each joint is actuated by a pair of agonist/antagonist muscles, each of which is individually activated by one neuron controller. The separate control of individual muscles mimics the distributed hierarchy of the neural control in real insect. Muscular Model. The typical Hill-type muscle model is introduced as the building block for the actuation system of the virtual insect 21,22 . To simulate a virtual muscle, we need to model the dynamics of activation and contraction.
The activation process converts the spiking trains from the neural system into the activation signal on the muscle. This is modelled as a first-order differential equation: where a t is the activation signal at time t, v t is the excitation from the motor neuron (Equation 2). Δt is the stepsize.
To simulate the contraction dynamics of a Hill-type muscle, we need to model two main features: the forcelength and force-velocity relationships. The force-length (FL) part models the function of both active and passive forces in terms of the muscle length, and the force-velocity (FV) part describes how the contraction force varies with respect to the muscle's contraction velocity.
where a t is the activation signal from Equation 6. The variables L ce , V ce are the length and contraction velocity of the muscle. To simplify the model, we set L ce = 1 in the equation of F v , which means that the change of muscle length is ignored when computing the force-velocity relationship. This assumption simplifies the model parameters as: are muscle-specific and determine the overall performance of a muscle. Existing works normally find the appropriate parameters using the collected data from vivo subjects 23 . However the experimental data are limited and may not be extended to different insects, or even different muscles on a single insect. Given the difference of the neural system between real and virtual insects, it is not of significance to follow the exact parameter values of real insects. Given the complicated model of nonlinear muscle, together with the parameters of neural system, we expect a high dimension of control space. It is challenging to reach an optimal solution for such optimization problem. In such circumstances, we propose the evolutionary algorithm to optimize the parameters of the controller and muscles simultaneously. Skeletal Model. During locomotion, an insect must solve essentially the same problem as a vertebrate although it has an external, rather than an internal skeleton. This is reflected in how the muscles are connected to the skeleton (Fig. 2). For the skeletal model, we follow the implementation from existing work 11 and model the virtual insect as a hierarchy of connected rigid body.
The main body is constructed as a set of three rigid bodies: head, thorax and abdomen, connected by two joints (head-thorax and thorax-abdomen). Each leg is actuated with three selected DOF: trunk-coxa, trochanter-femur and femur-tibia. Movement about the single joints and the resulting stepping patterns are generated by the activity of antagonistic muscle pairs. In the stick insect, the three major muscle pairs of a leg are the protractor and retractor coxae, the levator and depressor trochanteris, and the flexor and extensor tibiae. The protractor and retractor move the coxa, and thereby the leg, forward and backward. The levator and depressor move the femur up and down. The flexor flexes, and the extensor extends the tibia about the femur-tibia joint.
Each body segment is modelled as cylinder with its length listed in Table 1. The radius of the cylinder is set to uniform (0.2 mm) across all leg segments. The data are extracted from the anatomy of real ants. During the locomotion cycle, legs in stance mode are dynamically actuated and legs in swing mode are kinematically animated. This is supported by biologically evidences that each leg only occupy a small proportion of body weight (5%) ( Table 2), thus do not pose significant impact on overall body movements.

Optimization Algorithm
Previous paragraphs explain the hierarchy of our simulation framework. How to find the appropriate values for both the controller and muscle parameters is key to achieve stable performance of the locomotion. We here use the algorithm of Covariance Matrix Adaptation (CMA) to optimize the parameters. CMA is an evolutionary algorithm that is designed to solve nonlinear non-convex optimization problem. Existing works have applied this method to design the locomotion controller for virtual characters 11,12,28 . The optimization is initialized with the mean vector m and the standard deviation vector σ cma . For each new generation, λ offspring individuals are sampled with a normal distribution: i  The mean vector m represents the favourite solution from last generation. The covariance matrix C determines the shape of distribution ellipsoid. This covariance matrix is updated so that to increase the likelihood of successful individuals.
All individuals are evaluated and ranked by a fitness function f(x i ). μ individuals with the lowest fitness score are selected and named as the group of elites. The new mean of the next generation is updated from a weighted sum of these μ elites. The covariance matrix is updated similarly as a combination of covariance from previous generation. By iteratively repeating this process, the samples are expected to converge to the optimal solution.
In the model proposed in this work, the character has 3 DOFs (3 pairs of muscles) for each leg, and the dimensions of the neuron and muscle parameters are 3 and 8 respectively for each neuron-muscle pair. The symmetry between left and right legs first reduces the dimensions of motor and muscle parameters to half. Together with the synaptic weights in the central neural system, the total dimension of parameters is 806. This far exceeds the dimensions of parameter space in existing works 11,12,28 and the capability of the algorithm as declared by the author 29 . To optimize in the parameter space with such complexity, we resolve this challenge with a two-stage strategy with the algorithm of Covariance Matrix Adaptation (CMA). First, we optimize the input currents I s to the motor neurons and muscle parameters for the standard case of walking straight forward. Next, we optimize the offset currents from the central neurons, which are added to the standard case I s and generate the versatile motor skills including walking on curve trajectories and uneven terrains.
Objective Function for Stage I. The design of the objective function is critical in selecting the elites out of the whole population and optimizing the parameters via generations of evolution. The objective functions at Stage I are designed to consider three goals simultaneously.
The first is to move forward at a predetermined velocity ν * of body center (a vector pointing straight forward): where ν is the actual moving velocity of body center. The subscript k sums up all the simulation steps. The velocity is three dimensions (x, y, z) thus this objective function also penalizes the case when the body falls down or deviates from the forward direction.
The second is to match the foot trajectory p f with the captured trajectory ⁎ p f from the real insect: By minimizing this objective component, we expect to reproduce the similar gait pattern as the real insects. Ants typically demonstrate the double-tripod gait 11,30 . The intra-leg coordination between joints on the same leg is achieved by matching the trajectory of the particular leg, while the inter-leg coordination between different legs is enforced by the spatial-temporal information embedded in the trajectories of all legs on real ants. The use of motion trajectories allows the compact representation of the desired gait pattern, as part of the optimization goal.
The third is to minimize the energy consumption of the neural system, that is to use as little current input I as possible: This objective is biologically valid since the individuals use less energy for neural control gains advantages during the process of evolution. The final objective function is defined as: Objective Function for Stage II. The goal for the optimization in Stage II is to find the optimal values of synaptic weights in order to adapt to a variety of scenarios, specifically coping with curve trajectories and uneven terrains.
To accomplish this goal, we generate N s samples from 2D uniform distribution of two variables: and enforce an objective function focusing on the trajectories of body center and foot: By summing up the errors of all samples, we consider the fitness of the individual solution under general scenarios.

Results
Implementation. The code is written in C++ and runs on a standard PC with four CPU cores at 3.4 GHz.
The physics simulation is done with the Bullet 31 physics engine with a simulation step of 5 ms. Our single-thread implementation runs in real time during the online simulation, while the optimization process runs with the technique of multi-thread, costing around 10 hours for Stage I and 22 hours for Stage II on the aforementioned hardware. In comparison to the discretized look-up table 11 , our work provides a continuous mapping with comparable timecost. Our method also avoids manually selecting the value of the substep width which critically determines the quality and time performance of the look-up table method 11,32 . During the optimization, the simulation for each individual is run for 100 seconds (or 100 locomotion cycles) for Stage I, or 10 seconds (or 10 locomotion cycles) for Stage II, or is terminated in advance if the body falls down (the height of body center falls below a threshold). Figure 4 gives a screenshot of the locomotion sequence of a virtual ant in a variety of scenarios. . α starts from value of 0 and progressively increase to 1, with a step of 0.2, so that the initial value stays sufficiently close to the optimal value. Figure 5(b) shows the optimization progress for Stage II. The result shows that when the exploration range is small (for cases when α = 0.2 or 0.4), the optimization converges fast and smooth, while for other cases the optimization requires more iterations to converge.
Spiking pattern and muscle activation. Figure 6 plots the spiking pattern and muscle activation signals of the retractor neuron-muscle unit from the body-coxa joint on the left-sided and right-sided middle legs, when the character turns left. With the activation dynamics (Equation 6), the discrete spiking signal is converted into a continuous activation stimulus. The result shows that when the character turns left, the number of spike trains on the left-sided leg is greater than the one on the right-sided leg. This leads to a larger activation signal for the muscle actuation. The result also shows that the coordination between different legs is achieved by preferably adjusting the number of spike trains, instead of the amplitude of the spike.
Comparison of joint trajectory between the real and virtual ants. One of the objective functions is to minimize the difference between the trajectory of the real and virtual insects. We here compare the joint angles at three DOFs on the middle leg for both the real and simulated subjects (Fig. 7). We found that the overall shape of the trajectory is consistent between the real and simulated cases. However in the result of real cases, the joint angles decrease to a minor extent before increasing incrementally. This is not observed in the simulated case, in particular for the body-coxa joint. This mismatch may be caused by the posture adjustment of real ants, or possibly by the hardware limitation (image resolution and framerate) of data collection. In the simulation case, the body-coxa joints move at a strict forward pattern in a single motion cycle. Figure 8 plots the joint torques on the middle leg during the stance stage for both the simulated and real cases. The torque from the simulation cases are recorded as the torques applied on the joint constraints during the simulation step. The torque from the real cases are computed with the method of Virtual Model Control 11 . This method is based on Jacobian matrix, converting the ground force on the end-effector into the joint torques given the configuration of body hierarchy. The peak of the simulation torques arises in the first half of the stance stage for the body-coxa and femur-tibia joints, and in the middle of the stance stage for the coxa-femur leg. However, the muscle activation in the simulation case reaches its peak around the end of the stance stage, instead of the first half. This is possibly due to the fact that the joint rotation reaches its limit at the end of the stance stage, resulting in a small value of F L and F V (Equation 9) and thus a moderate value of joint torque. The torque profiles of real cases generally follow the shape of a parabolic curve and have the peak value around the middle of the stance stage. This is consistent with the ground force on the middle leg ( Fig. 9(b)), which presents a similar parabolic curve shape. Figure 9 plots the ground forces on the real ant in the Z (vertical up) direction during the stance stage. For the front leg, the force magnitude present small variations during the whole stance stage. For the middle leg, the force     magnitude increases initially, remains stable and decreases when the stance stage comes to the end. This creates a parabolic curve of the force profile. For the back leg, it presents a distinct peak during the first half of the stance stage, which is in consistent with the torque results in the simulation case.

Discussion
Advantages of spiking neuron model. Replicating the function of neural system in real organism is non-trivial, given the fact that the machinery underlying the observed behaviour is not fully understood 32 . Existing works have proposed different paradigms to model the behaviour neural system. For example, previous work uses a network of nonlinear oscillator to construct the Central Pattern Generator 11 . We believe the use of different models has its own pros and cons. The oscillator with inherent properties of limit cycle and attraction point could generate the basic pattern with no intervention from higher control system. The high-level neural system is encoded as a look-up table 11 , which is pre-computed and discretized. In comparison, the current work, which is based on spiking neuron models, provides a unified representation of neural system. The spiking pattern of the unit neurons can be used to seamlessly connect the procedures of sensory processing, high-level commands and low-level control strategy. Such a unified representation lays the foundation for processing sensory information via large-scale networks.
Advantages of Hill-type muscle. The use of nonlinear Hill-type muscle improves the naturalness of the synthetic motion on the examples of biped 28 , and research also show that it improves the 2D motion stability for insect simulation 33 . However the significance of non-linear muscle has not been fully evaluated in the example of a 3D insect model. We here compare the performance of two actuator choices: Proportional-derivative (PD) servo 11 and nonlinear Hill-type muscle, in terms of the naturalness of the synthetic motion (Fig. 7). The PD servo is defined as: are the rotational angles and velocity of the joint. k g , k d are the gain and damping coefficients. To quantitatively evaluate the difference between the ground-truth data and the simulated trajectories, we use the Euclidean distance between two vectors: T T  T  T T  T  T  T  [ , , ], are joint trajectories of rotation in one step cycle from the ground truth and simulated result respectively. T g is computed as the average of all collected trajectories of real ants.
The quantitative result of comparison (Table 3) shows that trajectories of our model presents less deviations from the ground-truth data than the model in 11 . This is mainly caused by the linear feature of the PD servo. We also notice that the peak values of the trajectories of the previous model 11 are generally higher than the ones of the proposed model in this work. The reason is possibly caused by the fact of over-shooting 34 . The problem of over-shooting refers to the fact that the PD servos normally require a large value of k g , thus causing the over-protraction/elevation/extension of the foot.
The reason is possibly caused by the fact of over-shooting 34 . The problem of over-shooting refers to the fact that the PD servos normally require a large value of k g , thus causing the over-protraction/elevation/extension of the foot.
The neuro-musculo-skeletal model we developed yields very promising results. The CMA optimization algorithm successfully helps to find optimal parameter configuration for the model. The optimization process converges quickly in around 2000 iteration. And the simulated gait pattern is highly consistent with the experimental data collected from real ants walking in the free mode.
An interesting and promising point is that the joint torques of the simulated and real cases also present high consistency. It will be a great help to biologist in studying of some special phenomenons of insects such as the incredible lifting ability of ants and impressive leaping ability of mantis. They can easily observe and analysis the stress conditions under various situations with the simulated model. And the study can stimulate the development of bionic design.
The model could also be beneficial in studying of higher level behavior of insects. With the model, we can reproduce the motion of an insect with only the recorded foot trajectories. Thus, we can use techniques of machine learning to extrapolate or predict the trajectories of insects under different circumstances and then reconstruct the 3d motion with our model. Our model is intended to be widely applicable in its ability to model insect locomotions. Here we only conduct a primary validation on ants. It will be interesting to apply this model to other insect species in future, potentially including other locomotion constraints.  Table 3. Quantitative comparison of the trajectory difference between the ground-truth data and the ones generated by our model and the previous model 11 . In conclusion, we have proposed a neuro-musculo-skeletal model of virtual insects. By incorporating some nature-tested mechanisms of real insects, the model is capable of reproducing the gait pattern observed from real insects on virtual ones. There are a couple of directions for future works. The current implementation uses a single neuron to simulate the neuron group activating the same muscle. Improving the model complexity to match the real insect is worth exploring. To extend this framework to these advanced motion skills is one of the future efforts. How to encode more sophisticated behaviours, such as collective transport 35 , as patterns of spiking trains is even more challenging.