Nonlinear neural control using feedback linearization explains task goals of central nervous system for trunk stability in sitting posture

Objective. Characterizing the task goals of the neural control system for achieving seated stability has been a fundamental challenge in human motor control research. This study aimed to experimentally identify the task goals of the neural control system for seated stability. Approach. Ten able-bodied young individuals participated in our experiments, which allowed us to measure their body motion and muscle activity during perturbed sitting. We used a nonlinear neuromechanical model of the seated human, along with a full-state feedback linearization approach and optimal control theory for identifying the neural control system and characterizing its task goals. Main results. We demonstrated that the neural feedback for trunk stability during seated posture uses angular position, velocity, acceleration, and jerk in a linearized space. The mean squared error between the predicted and measured motor commands was less than 0.6% among all trials and participants, with a median correlation coefficient r of more than 0.9. Our identified optimal neural control primarily used trunk angular acceleration and near-minimum muscle activation to achieve seated stability while keeping the deviations of the trunk angular position and acceleration sufficiently small. Significance. Our proposed approach to neural control system identification relied on a performance criterion (e.g. cost function) explaining what the functional goal is and subsequently, finds the control law that leads to the best performance. Therefore, instead of assuming what control schemes the neural control might utilize (e.g. proportional-integral-derivative control), optimal control allows the motor task and the neuromechanical model to dictate a control law that best describes the physiological process. This approach allows for a mechanistic understanding of the neuromuscular mechanisms involved in seated stability and for inferring the task goals used by the neural control system to achieve the targeted motor behavior. Such neural control characterization can contribute to the development of objective balance evaluation tools and of bio-inspired assistive neuromodulation technologies.


Introduction
Maintaining stability of the inherently unstable human trunk during sitting requires the complex interaction of several components of the sensorimotor system [1,2]. Individuals with moderate-to-severe sensorimotor impairments, e.g. up to two-thirds of individuals with spinal cord injury, exhibit reduced trunk stability and the inability to sit unassisted [3][4][5]. Consequently, they are at risk of injurious falls when exposed to external sitting disturbances, with nearly 70% of wheelchair users experiencing at least one fall each year [6]. Such falls are the leading cause of injury in this population [7] and often require hospital care [8]. In this light, previous work has demonstrated that an understanding of intact neuromuscular control can be utilized to identify and restore impaired balance in individuals with neuromuscular impairments [1,3,9]. Hence, the characterization of neuromuscular control can contribute to objective balance evaluation and the development of targeted rehabilitative interventions for improving impaired balance [10,11]. Furthermore, assistive technologies based on neuromodulation, such as functional electrical stimulation (FES), have shown potential for restoring seated stability by stimulating impaired trunk muscles via both open- [12] and closed-loop [7,[13][14][15][16][17] control strategies. However, designing bio-inspired controllers for FES systems is challenging [10,18,19], since such complex technologies must provide physiological actions similar to those of motor commands (e.g. muscle activation) produced by the intact neural control system [19]. Therefore, characterizing the neural control of nonimpaired seated stability is a prerequisite for bioinspired closed-loop neuromodulation technology [10].
To characterize neural control, it is essential to determine the task goals representing high-level functional goals of the central nervous system (CNS) that regulate seated stability. However, characterizing the task goals for a given motor behavior, such as seated stability, and how the CNS accomplishes these goals via neural control, has been a long-standing question in human motor control research [3]. The literature has proposed different task goals for maintaining upright standing posture based on the reduction of the body center of mass (COM) sway and pace of sway [2]. Nevertheless, determining the task goals for seated stability, their relative importance in the neural control strategy, and how the neural control regulates this motor task is still an open question. Answering this question requires characterizing the complex interrelation between the underlying neuromuscular stabilization mechanisms and a mechanistic understanding of the closed-loop control system for seated stability.
To characterize the closed-loop postural control system and the interrelation between neuromuscular mechanisms, previous studies have proposed neuromechanical models of human stability. The closedloop postural control system during sitting can be explained using an inverted pendulum with passive and active stabilization control mechanisms [2,10,[19][20][21]. Identifying the dynamics associated with the process of converting sensory information into motor commands through active feedback control by the CNS characterizes the neural control and its task goals. Previous studies have identified the neural control as well as the passive and active controls of standing and seated stability using linear closed-loop system identification techniques applied to body kinematics and muscle activation data recorded when the body was perturbed via external stimuli (e.g. moving support surface, external forces, or perturbed visual surround) [20,22]. Previous studies proposed that the neural control generates motor commands based on task goals of minimizing the angular position and velocity [1,2,23,24] as well as acceleration [1,2] with respect to the static upright posture. Many studies have provided non-parametric estimates of the neural dynamics in both standing and sitting postures using linear closed-loop system identification [1,2,10,25]. Other studies have provided parametric estimates of the neural dynamics by using linear controllers such as proportional-integral-derivative (PID) control [10,21,26,27], proportional-derivative (PD) control [22,23], and PD control with acceleration feedback [1,28] to model the neural control functioning for controlling postural stability.
Another approach has been proposed to use linear optimal control theory to characterize the task goals for motor behavior [2]. This theory characterizes task goals as the minimization of cost functions to obtain an optimal performance of the neural control system [29,30]. In contrast to classic models (e.g. PD or PID) for neural control, the optimal feedback control theory relies on a performance criterion allowing the motor task and neuromechanical model of the body to dictate a control scheme that best describes the actual physiological process and, thus, results in a control law that leads to the best performance [30]. Oftentimes, this approach assumes that the neural feedback is minimizing a cost function that penalizes large motor commands (e.g. muscle activations) and state variables (e.g. body angular kinematics) toward minimizing energy expenditure and guaranteeing stability, respectively. Parameterizing the cost function enables prioritizing the task goals of the neural control for motor behavior [2]. For a given neuromechanical model of the body, there is optimal neural feedback that minimizes this cost function. Hence, identifying the neural feedback control strategy for a given motor task (e.g. seated stability) would allow us to interpret the cost function and, therefore, explain the relative importance of the task goals [2].
Despite its potential, optimal control theory has not been utilized to characterize neural control and its task goals for regulating seated stability. Characterizing the neural control using PD or PID control schemes along with linear timeinvariant neuromechanical models ignore the timevarying nonlinear behavior of the underlying neuromechanical dynamics associated with physiological uncertainties in real-world conditions [31][32][33][34]. In addition, the use of control structures such as PID does not reflect optimal neural control strategies and physiological task goals for regulating seated stability. Recently, nonlinear control theory has been used to characterize the biomechanics of postural control in the sagittal plane [35,36]. Sultan et al [35] utilized a high-gain observer along with feedback linearization to simulate the behavior of the CNS for regulating standing stability in the presence of large perturbations and neural delays. Their results showed that the nonlinear control paradigm, along with physiological latencies, is important for understanding the behavior of the motor control process. Sultan et al [36] also investigated the use of a high-gain observer-based feedback linearization paradigm to simulate the CNS for regulating sit-to-stand transfer. They demonstrated that such control paradigm provides an accurate simulation of the underlying dynamic system and is reliable and effective for customizing human models. Therefore, utilizing a nonlinear neuromechanical model along with nonlinear control theory for identifying the neural control would lead to a better mechanistic understanding of the task goals for a given motor behavior (e.g. seated stability). Such an identification approach would enable researchers to overcome a fundamental challenge in human motor control, open new horizons in the mathematical identification of neuromechanical mechanisms, and contribute to objective balance evaluation, developing targeted rehabilitative interventions, and bio-inspired design of assistive technologies, all aimed to improve impaired balance.
This study aimed to address the above-mentioned issues by identifying the task goals of the intact neural control associated with seated stability. First, we used a nonlinear neuromechanical model of the seated human along with a full-state feedback linearization approach and optimal control theory for identifying the neural control. Second, we identified the parameters associated with a cost function that penalizes the motor commands and state variables to achieve seated stability via nonlinear neural feedback. Third, we characterized the task goals for seated stability by interpreting the identified cost function.

Participants
Ten able-bodied, young, and male individuals volunteered to participate in this study (age: 24 ± 4 years; weight: 76 ± 13 kg; height: 178 ± 8 cm; mean ± standard deviation). All volunteers reported no history of neuromusculoskeletal disorders or any impairments that may have affected their seated stability. Participants gave written informed consent to the experimental procedures, which were approved by the local research ethics board.

Data acquisition
We placed retro-reflective markers, bilaterally, on the anatomical landmarks [37] of the head, thoracic and lumbar segments of the trunk, arms, hands, pelvis, legs, and feet (figure 1). We also placed four markers on the four corners of the seat. We recorded the trajectory of the body and base-of-support during the experiments using a 12-camera motion capture system (Vicon Motion Systems Ltd, Oxford, UK) at a sampling rate of 100 Hz. The motion capture system (Vicon Motion Systems Ltd, Oxford, UK) was calibrated by the laboratory operator with the Vicon Calibration T-Wand, which featured five active redlight LEDs. The operator waved the wand around the capture volume where we intended to capture the 3D motion of retroreflective markers. The wand was finally placed at the corner of the platform-integrated force plate to set the origin of the global reference frame of the system.
We placed bipolar electromyography (EMG) electrodes (Bagnoli, DELSYS, Natick, MA, USA), bilaterally, over the muscle belly of the erector spinae (T9 level and L3 level), biceps femoris, rectus abdominis, external oblique, and rectus femoris (figure 1) [11]. The selected muscles are known for their major contribution to seated stability [38]. We used a selfadhesive reference electrode placed over the right iliac crest. We amplified EMG data with a muscle-and participant-dependent gain varying from 10 2 to 10 4 , sampled and digitized at 2 kHz.

Experimental procedure
We recorded the baseline muscle activity while the participants lay down in the supine position with their eyes closed for 60 s. We also recorded the maximum voluntary contraction (MVC) for each muscle while participants performed a series of exercises [11]. Each MVC was performed three times consecutively, with each time lasting 30 s. Participants were given a 30 s resting break between MVC exercises to minimize the effect of fatigue. A T-pose anatomical calibration trial was conducted (figure 1). The anatomical markers of the thoracic, lumbar, and pelvic segments were removed prior to the main trials. The anatomical calibration data were then later used to reconstruct the anatomical markers for the main trials.
We fastened a customized seat on a six-degreeof-freedom Stewart platform of a Computer-Assisted Rehabilitation Environment (CAREN; Motek Medical, Amsterdam, The Netherlands) before the session started. A consistent visual surround was provided for all participants by projecting two-dimensional, white grid lines on a black background on a 180-degree curved virtual-reality projection screen (figure 1). Participants sat on a customized seat with their arms crossed over the chest and their legs hanging vertically downward without foot support (figure 1). We asked the participants to keep their vision on the virtual-reality screen and maintain their seated stability while external disturbances were applied in the form of anteroposterior translation of the platform. First, we used a low-amplitude ramp-shaped perturbation suddenly applied to the seat. Subsequently, we used three 240 s white noise profiles as the perturbation. We designed the white noise profiles with a mean power spectral density of 4 cm 2 Hz −1 , filtered via zero-lag first-order high-pass and eight-order lowpass Butterworth filters with cut-off frequencies of 0.1 Hz and 5 Hz, respectively. Each profile was initiated with a 5 s hold followed by a 5 s increasing ramp and terminated with a 5 s decreasing ramp followed by a 5 s hold to eliminate any abrupt initiation and termination effects. Therefore, the middle 220 s of each trial were used for analyses. Between-trial rest breaks were given upon request by the participant.

Data pre-processing
We filtered the motion trajectories via a zero-lag fourth-order low-pass Butterworth filter with a cutoff frequency of 10 Hz. We estimated participantspecific body segment parameters including mass, moments of inertia, COM, and joint center of rotation positions by scaling cadaveric data of the upper body [39] and lower limbs [40] using each participant's weight and height [41,42]. Previous studies obtained segmental kinematics based on the marker trajectories [43,44]. We calculated the instantaneous position of the upper body COM using the weighted summation of the instantaneous COM position of the head and neck, thoracic, lumbar, arms, and pelvic segments. We used an inverted pendulum model to obtain the trunk sway angle as the motion of the upper body COM rotating about the L5-S1 joint.
We divided the EMG time series by gain, and we then demeaned, rectified, and filtered the scaled time series using a moving-average filter with a window size of 100 ms [45]. We down-sampled the EMG time series to 100 Hz and normalized them for the perturbation trials using a Min-Max Scaling approach based on the baseline and maximum MVC (highest values among the three recorded MVCs) of each muscle as the minimum and maximum values, respectively.

Nonlinear neuromechanical model
Previous studies have proposed a nonlinear neuromechanical model of seated stability [46]. The nonlinear neuromechanical dynamics of the seated stability was modeled, as illustrated in figure 2. A statespace representation of this model has four state variables, one input and one output. The first two state variables represent the angular position and velocity of the trunk modeled as an inverted pendulum controlled by gravitational, active, and passive joint moments. The other state variables represent critically-damped second-order transfer function associated with the muscle activation dynamics due to calcium release dynamics, muscle fiber conduction velocities, and the time delay associated with chemical reactions [32,47]. The state-space representation is as follows: (1) Figure 2. The nonlinear neuromechanical model of seated stability was validated in past research [46]. A state-space representation of this model has four state variables, and one input and one output. The first two state variables represent the angular position (x1) and velocity (x2) of the trunk COM controlled by gravitational ( where m is the total mass of the upper body; g is the gravitational acceleration; l is the length of the inverted pendulum, θ0 is the trunk angle deviation from the vertical axis during upright posture; K1, K2, and K3 are the stiffness coefficient, exponential elasticity, and resting elastic angle of the trunk, respectively [32,48]; B1 and B2 are the damping coefficient and exponential term, respectively [31]; and Fm is an arbitrary function (e.g. a polynomial or a neural network) modeling the maximal active moment in the case of full muscle activation based on the trunk angular position (x1) and velocity (x2). The other state variables (x3 and x4) represent a critically damped second-order transfer function [32,47,49] associated with the dynamic behavior of the muscle activation due to calcium release dynamics, muscle fiber conduction velocities, with the parameters ω0 as the natural frequency and β as a damping coefficient equal to unity as well as the time delay (T d ) associated with chemical reactions. Note that x3 is the normalized activation, whereas x4 is its rate of change used in the transfer function (equation (2)) depicted with a box in this figure. The total motor command (u) was calculated as the summation of the back muscles' motor commands (extensors) subtracted by the summation of the front muscles' motor commands (flexors): each muscle activation was normalized via Min-Max scaling using its maximum voluntary contraction (MVC) and baseline. External disturbances were modeled as a joint moment M dist .
In equations (1)-(4), f and g are the system functions; h is the measurement function;ẋ is the time derivative of the state-vector x; u is the total motor command; y is the system measurement; x 1 is the trunk sway angle; x 2 is the angular velocity; x 3 is the normalized activation; x 4 is the rate of change of activation; M g is the gravitational joint moment; M e is the passive elastic joint moment; M v is the passive viscous joint moment; M a is the active joint moment; M dist is the joint moment due to external disturbances; J is the moment of inertia; ω 0 is the natural frequency; β is a damping coefficient equal to unity; and T d is the input time delay.

Identification of passive-active controls
We used the nonlinear neuromechanical model to characterize the passive and active stabilization mechanisms (figure 2). We used the data of the sudden perturbation trial to estimate parameters of the passive stiffness (K 1 , K 2 , K 3 ) and damping (B 1 , B 2 ) by using a Nonlinear Least Squares algorithm by assuming that low-amplitude perturbations mostly elicit passive mechanisms rather than active postural control mechanisms [50]. Subsequently, we used an unscented Kalman filter (UKF) [51] on three 240-second perturbation trials to estimate the parameters of the active control in addition to the state variables for each participant using a dual-estimation scheme by forming an augmented state vector: where X is the augmented state vector; is the state-space representation of the augmented statespace model, and ∅ is the active control parameters including the natural frequency (ω 0 ) of the activation dynamics and the coefficients (ψ ) of the arbitrary function F m . It was assumed that ∅ was slowly varying compared to the process dynamics [32]. We have previously shown that such an identification approach, along with the nonlinear neuromechanical model, allows us to quantitatively characterize the ) to obtain full-state neural feedback. We used a full-state feedback linearization approach to linearize the nonlinear neural feedback and then employed optimal control theory to identify the actual neural control (shown in blue in both figures 2 and 3) and its task goals associated with seated stability. We assumed that the neural control acts as a linear quadratic regulator (LQR), and the neural feedback minimizes a quadratic cost function that penalizes the motor commands and state variables in the linearized space. We used a genetic algorithm to obtainQ andR that minimized the difference between the predicted motor command (û) and the actual motor command (u) measured by EMG.
active and passive postural control mechanisms with high accuracy [46].

Nonlinear neural feedback: full-state feedback linearization
The use of UKF, in the previous study, allowed us to observe the state vector (x) of the nonlinear system based on the measurement of the trunk sway angle (y) and the total motor command (u) [46]. In the present study, to characterize the task goals of the neural control, we assumed that: 1. the neural control system acts as an optimal linear quadratic regulator (LQR) that receives fullstate linearized feedback from the neuromechanical model, and 2. the neural feedback minimizes a quadratic cost function that penalizes the motor commands and poor stability performance expressed as the state variables in the linearized space.
Hence, identifying the neural control is to determine: 1. a diffeomorphism that transforms the state variables of the nonlinear neuromechanical model into a linear space where LQR can be applied. 2. The weights of the quadratic cost function associated with the LQR control; and 3. the gains of the LQR that minimizes the quadratic cost function.
This allows us to interpret the cost function, the optimal neural control, and therefore, explain the task goals used by the neural control for seated stability [2]. The process of identifying the neural control is depicted in figure 3, and the details are provided in the following sections. First, we prove that there exists a transformation that transforms the nonlinear neuromechanical model into a linear space where linear optimal control theory can be applied.

Proof of existence of a transformation
A nonlinear systemẋ = f (x) + g (x) u where f : D → ℜ n and g : D → ℜ n are sufficiently smooth in a domain D ⊂ ℜ n , is said to be full-state linearizable if a sufficiently smooth diffeomorphism T : D → ℜ n exists such that D z = T (D) contains the origin and a change of variables z = T (x) that transforms the systemẋ = f (x) + g (x) u into the formż = Az + Bv and u = φ (x) + W −1 (x) v, with (A, B) being controllable and W (x) nonsingular for all x ∈ D [52]. We obtained a diffeomorphism z = T (x) that transformed our nonlinear systemẋ = f (x) + g (x) u into the linear space with (A, B) controllable and W (x) nonsingular on a domain D ⊂ ℜ n , for all x ∈ D: where z and v are the transformed state vector and control command in linearized space, respectively; φ (x) and W (x) are the control mapping functions; A and B are linear system matrices in a canonical form, and n is the number of state variables (n = 4 in our study). g} is involutive in D 0 . Note, ad i f g represents the i-order Lie bracket of f and g [52]. We investigated these two conditions for our nonlinear neuromechanical model as follows: Based on our neuromechanical model, we can show:

Condition
Using our nonlinear neuromechanical model, we can obtain: g, ad f g, . . . , ad n−2 Based on equations (15)-(17), Condition (2) is also satisfied. Therefore, the existence of full-state feedback linearization is proved.

Linearized full-state feedback
We obtained the linearized full-state feedback as the input of the neural control by solving a set of partial differential equations for the diffeomorphism z = T (x). Combining equations (1)-(3) and (6)- (8) gives:ż Since equations (13) and (14) are equal, we can obtain: T (x) is not uniquely defined by this system of differential equations and any linear transformation M that leads toT (x) = M.T (x) will also satisfy this set of equations; however, withÃ = MAM −1 andB = MB. One possible solution of this set of equations can be obtained as follows [52]: where L i f h (x) is the i-order Lie derivative of h with respect to f. Interestingly, the linear state-space z represents the angular kinematics of the trunk with z 1 , z 2 , z 3 , and z 4 having the dimension of angular position [rad], velocity [rad s −1 ], acceleration [rad s −2 ], and jerk [rad s −3 ], respectively. Therefore, the input of the neural control would be the angular kinematics of the trunk.

Neural control identification
The process of identifying the neural control is depicted in figure 3. The UKF provided the state vector x based on the measured trunk angle y using the processed motion capture data and the measured total motor command using the EMG data for each participant and each trial. We used the transformation matrix T (x) (equations (17)- (20)) to obtain the state vector z in the linear space expressing the trunk's angular kinematics and fed it into the neural control defined as an LQR with a quadratic cost function. The neural control then output linearized control command v which was then converted into the motor commandû in the nonlinear space using control mapping functions, W (x) and φ (x), based on equations (7), (21) and (22).
The quadratic cost function penalized the motor commands (muscle activations), and linearized state variables as follows: Q is a diagonal four-by-four and R one-by-one positive definite matrices. To quantify the task goals associated with seated stability, we used a genetic algorithm to obtain Q and R that minimized the difference between predicted motor command (û) via the neural control model and the actual motor command (u) measured by EMG. The obtained Q and R were then used to calculate LQR gains according to equations (24) and (25) below. Algebraic Riccati equation (equation (2)) is solved to obtain S based on Q and R. The LQR gain vector (k) is obtained based on equation (25) (see figure 3). Note, equation (24) has multiple answers; however, the answer that makes the closed-loop system stable via k is selected. Also, inferring the obtained cost function via comparing the values of elements in Q and R allows us to characterize the task goals associated with seated stability,

Analysis of neural control model performance
We divided each 240 s trial into two sets to assess the performance of the identified neural control for each participant and trial: (1) an identification set, as the first 150 s of each trial, used to identify the neural control (Q, R); and (2) a test set, as the last 90 s of the trial, used to cross-validate the performance of the identified model. We used the Friedman test with a significance level of 0.05 to statistically compare the weights of the LQR quadratic cost function (Q, R) as well as the LQR gains among all participants and trials. We used Friedman test for statistical comparison since the normal distribution condition was not met for our data. We used the mean squared error (MSE) and correlation coefficient (r) between the actual measurements (u) and the output of the identified model (û) on the test set to quantify its performance. Note that the correlation coefficient (r) is unitless and can range from −1 to 1. A correlation coefficient of −1 indicates a perfect inverse correlation, with one time series increasing (decreasing) when the other time series decreases (increases). A correlation coefficient of 1 indicates a perfect direct correlation. As shown in the following section, in the present study, the correlation coefficients between u andû were positive, implying a direct correlation. Therefore, we multiplied the correlation coefficients with 100 to present the results in terms of the percentage. Finally, we investigated the robustness of the estimated neural control by deviating the neuromechanical parameters by 10% to see the effect of erroneous neuromechanical model parameters. We used the Friedman test with a significance level of 0.05 to find any significant differences between the predictions before and after deviating the parameters in terms of the accuracy and correlation.

Results
Figure 4(a) shows the normalized weights of the quadratic cost function associated with the LQR control for all participants and trials. The obtained cost function weights, Q and R, were normalized by the maximum absolute value of the associated state variables and motor command, respectively, for comparison. Q 1 , Q 2 , Q 3 , and Q 4 penalized poor angular position, velocity, acceleration, and jerk, respectively, while R penalized the motor command representing muscle activations. Friedman test showed that R was significantly larger than Q 2 and Q 4 (P < 0.0001). Among the state variables, we observed that Q 1 and Q 3 were significantly larger than Q 2 and Q 4 (P < 0.05). The LQR gains k 1 , k 2 , k 3 , and k 4 were normalized by the maximum absolute value of the linearized feedback: angular position, velocity, acceleration, and jerk, respectively. We observed that acceleration feedback had a significantly higher gain compared to other kinematics information ( figure 4(b)).
We observed that the MSE between the predicted motor command (û) and measured motor command (u) was less than 0.6% of the Min-Max scaled value among all trials and participants for both identification and test sets ( figure 5). The correlation coefficient betweenû and u was higher than 90% as median  figure 3, Q and R, were normalized by the maximum absolute value of the associated state variables and motor command, respectively, for comparison among weights. Q1, Q2, Q3, and Q4 penalized poor angular position, velocity, acceleration, and jerk, respectively, while R penalized the motor command representing muscle activations; (b) The linear quadratic regulator (LQR) gains k1, k2, k3, and k4 were normalized by the maximum absolute value of the linearized feedback: angular position, velocity, acceleration, and jerk, respectively. Statistical Friedman test was performed to identify significant differences at a significance level of 0.05.  among all trials and participants for both identification and test sets. The most and least accurate estimations of the motor command had a correlation coefficient of 99.3% and 82.02%, respectively ( figure 6(a)). In addition, figure 6(b) shows an example of the performance of the identified neural control and neuromechanical model compared to the measured inertial trunk moment for one participant and trial.
We investigated the effect of erroneous neuromechanical parameters on the performance of the identified neural control. To this end, we deviated the neuromechanical parameters by 10%. We observed no significant effect on the performance of the identified neural control in terms of the correlation coefficient and MSE (figure 7).

Discussion
This study aimed to characterize the task goals of the neural control system for seated stability by experimentally identifying this system using a nonlinear neural feedback model in non-disabled individuals. We obtained full-state feedback from a nonlinear neuromechanical model of seated stability via a UKF and then used a full-state feedback linearization along with optimal control theory to quantify the neural control and its task goals.

Task goals of neural control for seated stability
To infer the task goals, we assumed the neural control system acted as an optimal controller that receives full-state linearized feedback from the nonlinear neuromechanical model of seated stability. Full-state feedback linearization transformed the nonlinear state vector into linear space where the state vector included trunk COM angular position, velocity, acceleration, and jerk. Therefore, the quadratic cost function of the LQR penalized the motor commands (e.g. muscle activations) and COM angular kinematics representing sensory information.
Inferring the optimized weights of the cost function allowed us to determine the task goals for seated stability via the optimal control theory [2]. The obtained weights showed that the neural control had significantly higher penalties for the motor command, and the COM angular position and acceleration compared to the COM angular velocity and jerk ( figure 4(a)). This suggests that the identified neural feedback stabilizes the trunk while achieving nearminimum muscle activations as a task goal.
Kiemel et al [2] showed that the neural control adopted minimization of muscle activation for stabilizing upright stance, and the neural feedback did not substantially increase muscle activation to reduce deviations of the COM position. Significantly high penalty for the motor command shows that there is no functional benefit to further minimizing COM angular sway that would demand higher muscle activation when sitting is stable. In addition to what Kiemel et al [2] observed for standing stability, the high penalty associated with COM angular position in our study reflects that the neural control attempted to sufficiently minimize the deviation of the COM angular position from upright sitting posture while maintaining minimum muscle activation. Todorov [30] explained that such behavior is due to 'minimal intervention principle' meaning that there is no need to correct deviations away from the average behavior unless such deviations hinder the task performance. This is because correction is expensive due to the energy costs, i.e. muscle activation. Satisfying task goals of simultaneously minimizing both muscle activation and COM angular position seems to be conflicting in terms of energy expenditure, since minimizing COM position demands high muscle activation. However, the neural feedback law that maintains the deviations of COM sufficiently smaller than the size of the base of support can achieve seated stability while keeping the muscle activation to a minimum [2]. Such a stabilization strategy, while ensuring stable posture, avoids overcorrection of COM position and high muscle activation that could lead to rapid onset of muscle fatigue.
We also found that, in addition to the COM angular position, the neural control penalized the COM angular acceleration significantly more than COM angular velocity and jerk. This suggests that acceleration information plays a more important role than velocity information which contradicts previous studies that modeled the neural control via PD control with no acceleration feedback [22,23]. Particularly, one study by Jeka et al [24] suggested that velocity information is the primary input to stabilize posture during quiet standing, rather than position and acceleration. This different interpretation could stem from fundamental experimental and theoretical differences between the studies. First, their method was based on linear systems theory, modeling the COM angular trajectory by a linear dynamic model. The linearity assumption is valid for small COM motions seen in their experiments, where COM motion was considered equal to center of pressure motion. Such modeling ignores the inherent nonlinearity between the input (i.e. center of pressure as representative of the ankle torque) and output (i.e. COM motion) of the underlying system. Second, their model neglects the sensory time delay, and the extent of contributions of passive control to postural stability was neglected. As presented by previous studies [22,23,46], passive control primary depends on the angular position and velocity and has a significant contribution towards postural control. Therefore, their interpretation regarding the role of the velocity information for the neural feedback may have been affected by not separating the passive and active control mechanisms, since the passive control mechanism depends on the angular position and velocity. In contrast, the present study has distinguished the roles of the passive and active control mechanisms, allowing for isolating the contribution of the neural feedback towards postural stability.
The literature has highlighted that acceleration feedback provides reliable information on the trunk's inertia during stabilization task [53]. Furthermore, angular acceleration information provides sensory input about the sum of joint moments applied to the trunk. The literature has shown the COM acceleration during standing is highly correlated with the distance between the center of pressure and the vertical projection of the COM position [54]. This distance reflects the relationship between the controlling and controlled variables of postural control [54]. As a result, COM acceleration contains information on the neural control system's performance in correcting sitting or standing posture. Therefore, such information can be relevant and employed by the neural control to generate motor commands resulting in proper stabilizing active joint moments in response to external disturbances. Moreover, we observed significantly higher LQR gain associated with acceleration input compared to position and velocity ( figure 4(b)). This reflects a larger contribution of the trunk COM angular acceleration to optimal motor command and, thus, active joint moments, further highlighting the importance of acceleration feedback. One additional reason that COM angular acceleration had more contribution to the motor command compared to COM angular position could be due to minimizing the muscle activation without overcorrecting COM position that would demand higher muscle activation levels.
Overall, inferring the cost function of the optimal neural feedback revealed that the task goals of the neural control for seated stability are to achieve near-minimum muscle activation while keeping the deviations of the COM angular position and acceleration sufficiently small. To achieve seated stability, the neural control system predominantly uses COM angular acceleration information to generate motor commands resulting in activating relevant muscles and generating active joint moments against external disturbances while avoiding overcorrecting COM angular position and rapid onset of muscle fatigue.

Nonlinear neural feedback model and its validity
The inherently nonlinear dynamics of neuromuscular mechanisms involved in seated stability and their complex interrelations hinder our mechanistic understanding of how the CNS stabilizes the human trunk during sitting. Previous studies have assumed linear neuromechanical models of stability to quantify neural control by using linear controllers such as PID control [10,21,26,27], PD control [22,23], PD control with acceleration feedback [1,28], and LQR optimal control [2]. For approximation of the human body dynamics during sitting as a linear model, the COM motion must be small, which is not necessarily true under challenging conditions or in the presence of significant perturbations [2]. Therefore, neglecting the nonlinear behavior of the underlying neuromechanical system could lead to erroneous estimation of the neural control, which, consequently, affects our understanding of the task goals used by the neural control system to regulate stability. To address this challenge, we proposed a twostep approach. First, we used a validated nonlinear neuromechanical model of seated stability proposed by the literature [46]. Second, in the present study, we identified the nonlinear neural feedback control strategy in the CNS for seated stability by employing full-state feedback linearization technique along with a linear optimal control model, i.e. LQR.
In contrast to the literature, our identified neural feedback is not limited to a small range of motions and has the potential to be applicable for a wide range of COM motion states. This is due to the global transformation from the nonlinear state-space into a linear state-space where linear optimal control theory is globally valid. Nevertheless, further experimental investigation is needed to verify the validity of this model for larger motions. By modeling the nonlinearity, we showed that the identified neural feedback cannot be approximated using a PD control and linear neuromechanical models, previously reported in the literature [22,23]. We showed that the neural feedback depends not only on the COM angular position and velocity but also on the COM angular acceleration and jerk, while the COM position and acceleration information contributed significantly more to the motor command compared to its velocity and jerk. As such, using PD control with a linear neuromechanical model limits our ability to explain the task goals used by the neural control for seated stability. In contrast, the use of optimal control with quadratic cost function allowed us to infer the task goals associated with seated stability.
An optimal control scheme relies on a performance criterion (e.g. cost function) explaining what the goal is and subsequently, finds the control law that leads to the best performance [29,30]. Therefore, instead of assuming what control schemes the neural control might utilize (e.g. PD or PID), optimal feedback control allows the task and the neuromechanical model to dictate the control scheme that best describes the process of regulating seated stability [30]. This gives the optimal control a generality in practice [30]. Moreover, since the sensorimotor system consists of components with actions aimed to continuously enhance the behavioral performance (e.g. natural optimality), the optimal control can precisely describe the optimal nature of the sensorimotor function compared to PD or PID control schemes [29,30]. This gives the optimal control a theoretical advantage over classic alternatives (e.g. PD or PID) for describing the behavior of the neural control [30].
It should be noted that by using LQR control, we assumed the external perturbations used in the study do not alter the neural feedback law. This assumption was based on previous studies that suggested mechanical perturbations do not significantly change neural feedback law, in contrast to sensory perturbations, which may cause alteration in sensory reweighting and consequently, may change the neural feedback law [2,49].

Accuracy of the identified parameters of the neural feedback model
We examined the performance of the identified neural control for each participant and trial using MSE and correlation coefficient. The MSE between the predicted motor command and measured motor command (i.e. EMG recording) was less than 0.6% among all trials and participants, while the correlation coefficient between them was higher than 90% as median among all trials and participants for both identification and test sets (figure 5). The MSE values were within the range of accuracy of the measurement devices, suggesting that the identified neural feedback model was highly accurate. This accuracy can be seen in figure 6 for the most and least accurate estimations of the motor command as the output of the neural control system among all trials.
We also examined the robustness of the identified neural control by introducing error into the parameters of the identified neuromechanical model. We deviated neuromechanical parameters by 10%. The introduced error affects the feedback linearization and control mapping processes. We observed that the addition of these errors in neuromechanical parameters did not significantly change the performance of the identified neural control (figure 7). This shows that the identified neural controls were robust against potential variation in the neuromechanical properties (e.g. due to muscle fatigue). Note that we utilized UKF and dual-estimation scheme for identifying the properties of the nonlinear neuromechanical model as well as the state vector (figure 2 and equation (1)) [46]. This allowed for online correction of the neuromechanical parameters (e.g. using forgetting factors), which enabled tracking system variations due to muscle fatigue and physiological uncertainties [46].

Limitations
First, we only recruited able-bodied, young, and male participants highlighting the need for further investigation on other populations with different sex, age group, and/or neuromuscular impairments. Second, we scaled cadaveric anthropometric data by participant's height and weight to obtain individual-specific body-segment parameters, which may have introduced error into the neuromechanical parameters [41]. Third we used a one-segment model of the trunk, assuming intervertebral motions are negligible. Future studies should investigate the effect of multi-segment trunk models on seated stability [41,42]. Utilizing multi-segment models requires uncorrelated external perturbation profiles simultaneously applied to each segment level. Fourth, in the future studies, the external perturbation should also be applied in mediolateral directions to investigate seated stability in the frontal plane. Fifth, we used a quadratic cost function for the LQR, assumed as the neural control scheme. The type of the cost function could also be an optimization variable, among others. Nevertheless, the present study took a step toward understanding the neural control system's task goals associated with seated stability, which could provide the foundation for future work.

Conclusions
In this study, we experimentally identified a nonlinear model for neural feedback control of seated stability. By using a nonlinear neuromechanical model along with feedback linearization and optimal control theory, we inferred the task goals used by the neural control system to regulate seated stability. We showed that the neural feedback uses angular position, velocity, acceleration, and jerk in a linearized space. We observed that the neural control tries to achieve near-minimum muscle activation while keeping the deviations of the trunk COM angular position and acceleration sufficiently small. To achieve seated stability, the neural control uses COM angular acceleration information to activate relevant muscles to generate required active joint moments against external disturbances. Our proposed approach to identifying the neural feedback control allows for a mechanistic understanding of the neuromuscular mechanisms involved in seated stability while enabling inferring the task goals used by the neural control system to achieve the targeted motor behavior.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Code availability
The codes used in the current study are available from the corresponding author on reasonable request.
Steven Phan for their contributions to data collection and pre-processing.

Funding information
This work was financially supported by the Natural Sciences and Engineering Research Council of Canada (Grant: RGPIN-2016-04106) and the University of Alberta. Alireza Noamani was also supported by Alberta Innovates Graduate Student Scholarship.

Conflict of interest
None of the authors have potential competing interests to be disclosed.

Ethical statement
This study was approved by the research ethics board of the University of Alberta (Study ID: Pro00063998).