Profiles of movement speed and positional variability based on extended LQG for various noises

Stochastic optimal control has been studied to explain the characteristics of human upper-arm reaching movements. The optimal movement based on an extended linear quadratic Gaussian (LQG) demonstrated that control-dependent noise is the essential factor of the speed-accuracy trade-off in the point-to-point reaching movement. Furthermore, the extended LQG reproduced the profiles of movement speed and positional variability. However, the expected value and variance were computed based on the Monte Carlo method in these studies, which is not considered efficient. In this study, I obtained update equations to efficiently compute the expected value and variance based on the extended LQG. Using the update equations, I computed the profiles of simulated movement speed and positional variability for various amplitudes of noises in a point-to-point reaching movement. The profiles of movement speed were basically bell-shaped for the noises. The speed peak was changed by the control-dependent noise and state-dependent observation noise. The positional variability changed for various noises, and the period during which the variability changed differed with the noise type. Efficient computation in stochastic optimal control based on the extended LQG would contribute to the elucidation of motor control under various noises.

www.nature.com/scientificreports/ and feedback controller based on the extended LQG 11 . The controller reproduced the human-adapted complex reaching trajectory. However, in the studies mentioned above about the extended LQG, the expected value and variance were computed based on the Monte Carlo method. The Monte Carlo method computes the expected value by repeating many trials and requires much computation time. It is important to efficiently obtain the expected value and variance to investigate whether the optimal feedback control can reproduce the characteristics of human movement in environments related to the effects of various noises. In this study, I obtained update equations to efficiently compute the expected value and variance based on the extended LQG. Using the update equations, I computed the profiles of simulated movement speed and positional variability for various amplitudes of noises in a point-to-point reaching movement. In the computations of simulated movements, control-dependent noise is always considered, while the other noise is considered from the state-dependent observation noise, state noise, state-independent observation noise, and internal estimation noise. Figure 1 represents the simulated movements based on the extended LQG for the coefficient σ c of the controldependent noise. Figure 1a shows the speed profiles. The profiles were almost bell-shaped. The speed peak was about half the movement duration when the noise was absent (i.e., σ c = 0 ). Although the peak appeared earlier when the noise was slight ( σ c = 0.2, 0.5 ), and when the noise was larger ( σ c = 1.0, 1.5, 3.0 ), the peak appeared later. When the noise was further larger ( σ c = 3.0 ), the velocity at the movement's end was far from zero. Figure 1b shows the profiles of positional variability. When the noise was absent ( σ c = 0 ), the standard deviation of the position was zero (i.e., the line was identical with the X-axis). When the noise was larger, the variability peak appeared later and the positional variability was larger over the movement. When the noise was even larger ( σ c = 3.0 ), the peak disappeared (i.e., the variability increased monotonically). The optimal movement when σ c was 1.5 was the result for the basic parameter set. In the basic parameters, only the control-dependent noise was considered. The profiles of movement speed and positional variability resembled the characteristics reported by past studies [2][3][4] . Figure 2 represents the simulated movements based on the extended LQG for the coefficient σ d of the statedependent observation noise. Figure 2a shows the speed profiles. All profiles were almost bell-shaped. As the noise became larger, the speed peak appeared later. Figure 2b shows the profiles of positional variability. Unlike the case of control-dependent noise, no significant change occurred up to about half of the movement, although the variability was larger in the latter half of the movement. Figure 3 represents the simulated movements based on the extended LQG for the coefficient r of the control cost. Figure 3a shows the speed profiles. As the control cost was larger, the speed peak appeared later. Figure 3b shows the profiles of positional variability. Unlike the case of state-dependent observation noise, there was a large change up to about half of the movement but little change after half the movement. In other words, as the ratio of state cost increased, the variability around the middle of the movement became larger. Figure 4 represents the simulated movements based on the extended LQG for the coefficient σ ξ of the state noise. Figure 4a shows the speed profiles. The profiles were mostly bell-shaped and largely remained unchanged by the state noise. This characteristic was also confirmed for state-independent observation noise ( Fig. 5a) and internal estimation noise (Fig. 6a). The scaling factor of the state-dependent observation noise in these results was zero ( σ d = 0 ). Although the resulting graph is not shown, the velocity peak appears slightly later when the scaling factor is not zero. Figure 4b shows the positional variability. As the state noise became larger, the variability increased over the movement. On the other hand, the variability in the state-independent observation noise and internal estimation noise did not change almost from the beginning to the middle of the movement (Figs. 5b, 6b).  Stochastic optimal control has been studied to explain the characteristics of human upper-arm reaching movement. In particular, the extended LQG framework has been verified for movement tasks, including various disturbances 5,8,9,11 . The disturbances can be treated as noise in the extended LQG. However, in the previous studies, the expected value and variance of the optimal movement were computed based on the Monte Carlo method. The Monte Carlo method computes the expected value by repeating many trials and requires much computation time. In this study, the updated equations provided an efficient computation of expected value and variance. Efficient computation in stochastic optimal control based on the extended LQG would contribute to the elucidation of motor control under various disturbances.

Methods
I computed the movement speed and positional variability based on the extended LQG 8 for various amplitudes of noises. Instead of using Monte Carlo method, I obtained update equations to compute the mean and variance of state variables.
Framework of extended LQG. The linear movement dynamic system in discrete time t is given as: www.nature.com/scientificreports/ where A and B are the dynamic system matrices, x is the state variables, and u is the control signal. State noise ξ is state-independent noise which is Gaussian with mean 0 and covariance � ξ ≥ 0 . The noise ε related to controldependent noise is Gaussian with mean 0 and covariance � ε = I . C is the scaling matrix for control-dependent noise. The matrix C is BF, where F is the scaling factor. It is already known that the initial state is mean x 1 and covariance 1 .
The feedback system is as given below: where y is the feedback signal and H is the observation matrix. State-independent observation noise ω is Gaussian with mean 0 and covariance � ω ≥ 0 . The noise ǫ related to state-dependent observation noise is Gaussian with mean 0 and covariance � ǫ = I . D is the scaling matrix for state-dependent observation noise.
The state estimate x is calculated as given below: where K is the filter gain matrix and internal estimation noise η is Gaussian noise with mean 0 and covariance � η ≥ 0.
Cost per step is calculated as given below: where Q is the matrix of the state cost and R is the matrix of the control cost. The optimal control u is recursively computed in the opposite direction over a period of time, as follows: where L is the control gain matrix. It is initialized as follows: S x n = Q n , S e n = 0 , and s n = 0 . Total expected cost is calculated as follows: The optimal filter is recursively computed forward in time as follows: It is initialized as follows: e 1 = 1 , x 1 = x 1 x T 1 , and xe 1 = 0.
Mean and variance of state variables. The update equations used to compute the mean and variance of state variables were obtained from Eq. (1) and the following equations.
(1) www.nature.com/scientificreports/ where the estimation error e t is x t − x t . Equations (10) and (11) are derived from Eqs. (2) and (3), and Eqs. (1) and (10), respectively. The mean of state variables can be computed sequentially by applying the following update equations.
The variance of state variables can be computed sequentially by using the following update equations: Application to the reaching movement. The simulated movement task was a single joint reaching movement similar to the study 8 . The movement was replaced with a translational point-to-point reaching movement for simplicity. The movement duration t end was set to 0.5 s. The starting position was 0 m, while the target position p * was 0.2 m. where t is the time step, m is the point mass, and τ 1 and τ 2 are the time constants (Table 1). The state x of the system is [p(t),ṗ(t), f (t),ḟ (t), p * ] T , where p is the position and f is the force. The initial state was x 1 = [0, 0, 0, 0, p * ] T and 1 = 0.
The matrices of the feedback system are depicted below: where σ d is the scaling factor. Table 2 shows the parameters of the feedback system. The state noise, state-independent observation noise, and internal estimation noise have covariances � ξ = σ 2 ξ I , � ω = σ ω diag[ω p , ω v , ω f ] The simulated movements were obtained by changing each value of σ c , σ ξ , σ ω , σ d , σ η , and r in the basic parameter set ( Table 3). The optimization was completed when the absolute value of the relative change in the total expected cost became less than the convergence tolerance ( 1.0 × 10 −15 ).