Andrii Shachykov
email: andrii.shachykov@loria.fr

Patrick Henaff
email: patrick.henaff@loria.fr

Anton Popov

Alexander Shulyak
email: shulyak.alex@mail.ru

Neuro-musculoskeletal Simulator of Human Rhythmic Movements

Keywords: central pattern generator, musculoskeletal model, Parkinson's disease I

In this article, a new neuro-musculoskeletal simulator of human rhythmic movements is presented. It is based on a specific model of central pattern generator able to generate variable rhythmic signals for controlling a biological joint by changing its intrinsic parameters from upper controller. The CPG consists of three layers and four types of neurons and controls human leg hip joint. Its output signals generated by motoneurons are applied to excitation inputs of two modelled muscles of the human hip. The variation CPG speed showed short transition process with stabilization on two leg motions in swing phase. Phase diagram of hip joint angle showed stable movement as stabilization on single trajectory. Further work is aimed on development of stable human gait by applying CPG for each leg muscle.

INTRODUCTION

Parkinson's disease (PD) is a neurodegenerative disorder marked by tremor, muscular rigidity, and slow, imprecise movement mainly affecting individuals over age of 60 [START_REF] Hirsch | The Incidence of Parkinson's Disease: A Systematic Review and Meta-Analysis[END_REF]. This disorder results primarily from the death of dopaminergic (DA) neurons in the substantia nigra pars compacta [START_REF] Dauer | Parkinson's Disease: Mechanisms and Models[END_REF]. Loss of DA leads to increase of required effort to move. Among other influence, PD affects control of movements, especially upper and lower limbs, causes tremor and rigidity [START_REF] Dauer | Parkinson's Disease: Mechanisms and Models[END_REF]. Some authors propose computational models of altered walking gait [START_REF] Muralidharan | A computational model of altered gait patterns in parkinson's disease patients negotiating narrow doorways[END_REF] or altered grip gesture [START_REF] Gupta | Computational model of precision grip in Parkinson's disease: a utility based approach[END_REF]. However, these works does not model or simulate lower neural structures that receive signals from basal ganglia (BG). These structures located in spinal cord are called central pattern generators (CPG) and are involved in rhythmic movements and control muscle activities [START_REF] Brown | On the nature of the fundamental activity of the nervous centres; together with an analysis of the conditioning of rhythmic activity in progression, and a theory of the evolution of function in the nervous system[END_REF].

A global model should consist of three major parts or levels:

 first level: a computational model of basal ganglia whose output mediates the second level;

 second level: a model of spinal structures composed of several CPGs projecting to the muscles through motoneurons;

 third level: a musculoskeletal simulation of human lower limbs that will execute locomotor movements including physical effects of articular multi-bodies systems interacting with the ground.

Aiming to simulate impact of PD disorders on human walking gait, this paper presents a neuro-musculoskeletal simulator for the second and third level of this global model. This paper is structured as follows. Section II briefly tells about state of the art in fields of modelling subcortical structures, CPGs, and musculo-skeleton system. Section III describes used models of CPG and musculoskeletal system. Section IV presents simulation results obtained for the human hip articulation. Finally, Section V concludes this work, with some remarks and future work directions.

II. STATE OF THE ART

A. Central pattern generator

Central pattern generator is a set of inter and motoneurons located in the spinal cord. A century of evidence has led to the conclusion that rhythmic locomotion activities are largely controlled by this network of spinal neurons [START_REF] Brown | On the nature of the fundamental activity of the nervous centres; together with an analysis of the conditioning of rhythmic activity in progression, and a theory of the evolution of function in the nervous system[END_REF][START_REF] Guertin | Central Pattern Generator for Locomotion: Anatomical, Physiological, and Pathophysiological Considerations[END_REF]. They can generate patterns automatically by themselves, without descending signals from the higher centers. However, sensory feedback controls the shape of generated patterns and descending signals contribute to the synchronization of inside and between CPGs [START_REF] Amrollah | On the Role of Sensory Feedbacks in Rowat-Selverston CPG to Improve Robot Legged Locomotion[END_REF].

Current usage of CPG in research, modelling, and application vary from investigation of purpose and regimes of specific groups of neurons in different segments of spinal cord [START_REF] Kiehn | Locomotion: Circuits and Physiology[END_REF][START_REF] Marder | Central pattern generators and the control of rhythmic movements[END_REF][START_REF] Guertin | The mammalian central pattern generator for locomotion[END_REF] through modelling neural networks resembling CPG [START_REF] Rowat | Learning algorithms for oscillatory networks with gap junctions and membrane currents[END_REF][START_REF] Rybak | Modelling spinal circuitry involved in locomotor pattern generation: insights from deletions during fictive locomotion[END_REF][START_REF] Ijspeert | Central pattern generators for locomotion control in animals and robots: A review[END_REF] to synthesis of control units in robotics that have the same behavior patterns of biological CPGs [START_REF] Nassour | Multi-layered multi-pattern CPG for adaptive locomotion of humanoid robots[END_REF][START_REF] Manoonpong | Adaptive, Fast Walking in a Biped Robot under Neuronal Control and Learning[END_REF][START_REF] Righetti | Programmable central pattern generators: an application to biped locomotion control[END_REF].

Model of CPG in this work is our own and was used before [START_REF] Nassour | Multi-layered multi-pattern CPG for adaptive locomotion of humanoid robots[END_REF]. It is based on work of Rybak et al. for a two-level CPG that separates the timing and activation of the locomotion cycle [START_REF] Rybak | Modelling spinal circuitry involved in locomotor pattern generation: insights from deletions during fictive locomotion[END_REF].

B. Subcortical structures

The models of subcortical structures (sometimes including motor cortex) could be a one-compartment conductance-based Hodgkin-Huxley [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] 3-neuron system to look for the origin of parkinsonian tremor [START_REF] Dovzhenok | On the Origin of Tremor in Parkinson's Disease[END_REF] or more complex firing rate model [START_REF] Pavlides | Computational Models Describing Possible Mechanisms for Generation of Excessive Beta Oscillations in Parkinson's Disease[END_REF]. Those models result in firing rates of neurons that, for the current purpose, need to be transformed into control commands for CPGs. Other models are more applicable and simulate handwriting [START_REF] Contreras-Vidal | A neural model of basal ganglia-thalamocortical relations in normal and parkinsonian movement[END_REF], object lifting [START_REF] Gupta | Computational model of precision grip in Parkinson's disease: a utility based approach[END_REF], and walking [START_REF] Muralidharan | A computational model of altered gait patterns in parkinson's disease patients negotiating narrow doorways[END_REF].

C. Musculo-skeleton simulation

To simulate a musculo-skeleton system, one could start from a cellular level [START_REF] Röhrle | A Physiologically Based, Multi-Scale Model of Skeletal Muscle Structure and Function[END_REF], but that would be an overdoing as an aim hasn't been for the most realistic simulation of human lower limbs; any simulator that uses classic mechanical Hill muscle model [START_REF] Hill | The Heat of Shortening and the Dynamic Constants of Muscle[END_REF] will do. Such model, based on work of Thelen [START_REF] Thelen | Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults[END_REF] is implemented in OpenSim (http://opensim.stanford.edu). But for initial research, a simpler and faster model was used [START_REF] Ackermann | Optimality principles for model-based prediction of human gait[END_REF][START_REF]Musculoskeletal Model for Simulation of Walking -Human Motion and Control Laboratory[END_REF], that is described further.

III. MATERIALS & METHODS

A. Model of one joint CPG

The CPG model used in this paper was initially proposed in previous work where it was used as a controller for the walk of a humanoid robot [START_REF] Nassour | Multi-layered multi-pattern CPG for adaptive locomotion of humanoid robots[END_REF]. Although this model is rather mathematical, it is supported by two neurophysiological studies and combines their propositions in multi-layered multipattern CPG model. This model uses signal from high-level controller (e.g. mesencephalic locomotor region (MLR)) that influences the frequency of generated patterns, but this signal can be generated by model of basal ganglia. The CPG architecture is composed of three layers (Fig. 1) and four types of neurons. In this work, sensory neurons are connected only to motoneurons.

1) Rhythmic layer (RG):

RG neurons are based on a neural model proposed by Rowat and Selverston [START_REF] Rowat | Oscillatory Mechanisms in Pairs of Neurons Connected with Fast Inhibitory Synapses[END_REF] with self-rhythmic generation ability, able to oscillate according to two parameters of membrane conductivities for fast and slow currents. With different values of the cell parameters, RG neurons can generate different patterns: quiescence, almost an oscillator, endogenous oscillator, plateau, depolarization, and hyperpolarization.

The model cell is represented by two differential equations, for the cell membrane potential and the lumped slow current:

     ,, m f inj dV fast V q i dt           , m dq q q V dt        
with τm<τs, the time constants of the neuron membrane and slow current activation.

  

, tanh , where V is the membrane potential, q the lumped slow current, τm the membrane time constant, τs the slow current time constant for activation, iinj the injected current. fast(V, σf) is an idealized current-voltage curve for the lumped fast current, σf is dimensionless shape parameter for the currentvoltage curve, q∞(V) is the steady-state value of the lumped slow current, Es is a reversal potential, σs is the potassium conductance normalized to leak conductance. The most influential parameters are σs, σf, and iinj (more details in [START_REF] Nassour | Multi-layered multi-pattern CPG for adaptive locomotion of humanoid robots[END_REF]).

f ff f fast V V A V A                     , ss q V V E       

2) Pattern-formation layer (PF):

PF neurons shape patterns generated by RG layer and chose which rhythm (flexion/extension) will dominate for a joint. They are also capable of rhythm deletion without resetting the phase of RG layer. As result, it is possible to deactivate motoneurons while RG continue to oscillate. As an activation function for PF neurons, a sigmoid function is used, so these neurons are managed by two parameters of amplitude and saturation.

     1 , 1 MLR MLR i I PF e              1 , 1 n rg pf i j j j w RG w S I n          
where PFi is the activation value of the i th pattern-formation neuron, α is a value that denotes the slope of the sigmoid function, θ is the center point of the curve that denotes the threshold of the neuron, I is the average input to patternformation neurons, wrg→pf is the weight of the synaptic connection between RG neurons and PF neurons, RGi is the activation of the i th rhythm generator neuron, Sj is the activation of the proprioception or exteroception neuron, and wj is the weight between this neuron and the pattern-formation neuron. αMLR is a single value that represents the descending control from the high-level controller. θMLR is the modulation of the threshold by the high-level controller that drives the rhythm domination (extension/flexion).

3) Motor layer (MN):

MN layer neurons locally control the muscles with input from previous PF layer and proprioceptive SN. Proprioceptive and exteroceptive SN measure angular position of joint or contact force of foot and inhibit the corresponding motoneuron thus implementing articular reflex. MN and SN also use sigmoid activation function.

   1 , 1 i I MN e        1 , 1 n pf mn i j j j w PF w S I n          
The model also includes sensory neurons (SN) that shape the activity of the CPG neurons:

   1 , 1 i e         
where 𝜙 represents the joint angle.

A. Musculoskeletal model

In this paper the framework "GAIT2DE" developed by Ton van den Bogert, Orchard Kinetics LLC [START_REF] Ackermann | Optimality principles for model-based prediction of human gait[END_REF][START_REF]Musculoskeletal Model for Simulation of Walking -Human Motion and Control Laboratory[END_REF] is used. Gait2de is a dynamic model that simulates muscle activities and their actions on skeleton to produce movements in the sagittal plane. It has nine kinematic degrees of freedom, seven body segments, sixteen muscles [START_REF] Hill | The Heat of Shortening and the Dynamic Constants of Muscle[END_REF], and its dynamics and outputs are differentiable with respect to all inputs. This model is implemented as Matlab MEX function and it takes ~0.03 ms to compute. Body segments of global model are trunk, and thigh, shank, foot in each leg. Each of them has the following parameters: mass, length, center of mass, moment of inertia of a male with body mass 75 kg and body height 1.8 m. This simulator uses an implementation of well-known Hill muscle model [START_REF] Hill | The Heat of Shortening and the Dynamic Constants of Muscle[END_REF] that consists of parallel elastic element (PEE), damping element (DE), series elastic element (SEE), and contractile element (CE) (Fig. 2). Only CE is controlled by neural excitation u(t).

           1 () , act deact u t u da t u t a dt T t T t          
where a(t) is an active state of a muscle, Tact and Tdeact are time constants for activation and deactivation. More in [START_REF]Musculoskeletal Model for Simulation of Walking -Human Motion and Control Laboratory[END_REF].

The force in CE, or muscle fibers, depends on its active state, length, and lengthening velocity. The linear DE ensures that muscle model is solvable when the force in the CE is at or outside the asymptotic values of the force-velocity relationship. PEE and SEE represent passive properties of muscle fibers and surrounding tissue.

IV. SIMULATION RESULTS

As a preliminary step of our future work, we propose to simulate the neuro-musculoskeletal system of the right human hip (Fig. 3) that is the primary joint involved in the human walk. The two muscles mainly involved in hip movements are Iliopsoas and Glutei, and their parameters are given in Table 1. Since this work tells about control of only one joint, trunk was fixed at coordinates (0; 1) and 0 rad angle so it will not fall, interact with ground, or affect model with its swinging mass. To simulate the neuro-musculoskeletal system of the human hip, we had to modify the Gait2de simulator in order to send the signals generated by the CPG to the muscle model.

A. Connecting CPG to the muscles

Overall, full scheme of simulator (Fig. 4) consists of MLR projection to CPG described on Fig. 1, particularly as sinusoid iinj and σs, αMLR, and θMLR parameters. iinj enables RG desynchronization and could provide its rate for RG to adopt [START_REF] Nassour | Multi-layered multi-pattern CPG for adaptive locomotion of humanoid robots[END_REF]; σs controls RG frequency, it is actually a cell parameter, but it can be affected from upper controller; αMLR and θMLR control PF neurons, their coupling to MN and balance between flexion and extension.

The input to Gait2de are neural excitations (uF and uE) for each muscle, along with initial state of the model and optional external forces and moments applied to body parts and joints. MN provide neural excitations for right leg hip muscles (Iliopsoas and Glutei (Table 1)) that control hip joint. [START_REF]Musculoskeletal Model for Simulation of Walking -Human Motion and Control Laboratory[END_REF]) Iliopsoas is the hip flexion muscle that turns thigh forward, Glutei are three hip extension muscles, that are used as one here and turn thigh backwards. Additionally, output from MN is limited to [0; 1] as model does not apply such itself.

Muscle muscle mechanics muscle-skeleton coupling

To close the control loop, CPG contains SN for articular reflex [START_REF] Nassour | Multi-layered multi-pattern CPG for adaptive locomotion of humanoid robots[END_REF] that transform angle of corresponding joint into inhibitory influence on each motoneuron.

B. Constant speed of leg swing

First of all, an important feature of CPG is that it is able to produce rhythmic output without any input (Fig. 5). Tiny pulse could be required to desynchronize in case of same initial values both half-centers.

Next, outputs of RG half-centers are connected to PF neurons (Fig. 6), which in this case simply transform input value range to [0; 1], as flexion/extension domination and rhythm deletion through changing αMLR and θMLR aren't applied.

PF outputs are connected to MN, whose parameters are set to make joint movement symmetrical and softly limited to ±0.5 rad by SN. Flexor MN peaks (connected to Iliopsoas) are bigger than those of extensor because Iliopsoas muscle is weaker than Glutei (Table 1).

MN outputs excitate muscles which rotate hip joint. Muscle's variable length of CE and force are shown on the bottom of Fig. 6. SN react to maximum values of joint angle and inhibit MN thus preventing their turning above desired values.

This control result in stable joint rotation (Fig. 67). Change of speed on minimum and maximum angles happens because of physical joint limitation on minimum angle and interacting with free-swinging knee on maximum angle.

C. Changing the speed of leg swing

Speed of RG is controlled by neural parameter σs which always should be a positive value. Changing σs online is equivalent to reproducing neural plasticity effects. Indeed, σs can be modified by a learning or adaptive law modelling a homeostatic effect. After σs has changed, RG need some time to stabilize (Fig. 89). After speed change, joint produces one swing with higher amplitude and then stabilizes on second trajectory. This high-amplitude swing is produced by instantly changed σs parameter. It shouldn't be changed like that in natural control systems and is made in sake of studying the CPG controller. Other way is integrating several CPGs for each antagonistic pair of muscles for both legs thus achieving stable walk, driven only by CPG. However, this requires implementing some external means of achieving equilibrium. It is essential for simulating impact of PD disorders on human walking gait.

In addition, the last way may be in upgrading the physical simulator and musculoskeletal model to be more physically and biologically advanced.

Fig. 1 .

 1 Fig. 1. Model of one joint CPG controller with three layers: rhythm generator (RG), pattern-formation (PF), and motoneuron (MN) layers (extracted from [14])

Fig. 2 .

 2 Fig. 2. Arrangement of the elements in the muscle model. The CE is the only element that is controllable (extracted from [29] and modified). LCE is length of contractile element and LM is full muscle length

FmaxFig. 3 .Fig. 4 .

 34 Fig. 3. Scheme of Gait2de model used in this work

Fig. 5 .Fig. 6 .

 56 Fig. 5. Rhythm generator with no input. F for flexion RG half-center, E for extension

Fig. 7 .

 7 Fig. 7. Phase diagram of hip joint control from Fig. 6

Fig. 8 .

 8 Fig. 8. RG speed changing at 30 seconds and its effect on the system

Fig. 9 .

 9 Fig. 9. Phase diagram of hip joint control from Fig. 8. Solid line for initial higher speed before 30 seconds, dashed line for lower speed after 30 seconds V. CONCLUSION This paper presents a new neuro-musculoskeletal simulator of human rhythmic movements, aiming to simulate impact of PD disorders on human walking gait. A manner of controlling musculoskeletal joint with a model of central pattern generator is described. The simulator is based on connecting model of central pattern generator to a simulator of human musculoskeletal system. It is able to generate rhythmic movements for controlling a biological joint and variate their parameters from upper controller. The CPG consists of three layers and four types of neurons: rhythm-generation, patternformation, and motoneuron layers plus sensory neurons for articular reflex. CPG controls two muscles, Iliopsoas and Glutei, of human leg hip joint. Phase diagram of joint angle showed stable movement as stabilization on single trajectory. The variation of CPG rhythmic activity showed short transition process with stabilization on two trajectories.Further work may be aimed on several ways. One way is development of more complex CPG model to study more afferent feedback sensors and connections.

TABLE I .

 I PARAMETERS OF MUSCLES USED IN THIS EXPERIMENT (EXTRACTED FROM

ACKNOWLEDGMENT

This work is supported by Erasmus+ 2016-2017 and French embassy of Ukraine.


published or not. The documents may come

L'archive ouverte