Quantitative evaluation of posture control in rats with inferior olive lesions

Impairment of inferior olivary neurons (IONs) affects whole-body movements and results in abnormal gait and posture. Because IONs are activated by unpredicted motion rather than regular body movements, the postural dysfunction caused by ION lesions is expected to involve factors other than simple loss of feedback control. In this study, we measured the postural movements of rats with pharmacological ION lesions (IO rats) trained to stand on their hindlimbs. The coordination of body segments as well as the distribution and frequency characteristics of center of mass (COM) motion were analyzed. We determined that the lesion altered the peak properties of the power spectrum density of the COM, whereas changes in coordination and COM distribution were minor. To investigate how the observed properties reflected changes in the control system, we constructed a mathematical model of the standing rats and quantitatively identified the control system. We found an increase in linear proportional control and a decrease in differential and nonlinear control in IO rats compared with intact rats. The dystonia-like changes in body stiffness explain the nature of the linear proportional and differential control, and a disorder in the internal model is one possible cause of the decrease in nonlinear control.


Supplementary Methods
Brain sections of rats After all of the experiments, the rats were deeply anesthetized, their chest opened, and physiological saline and 4% paraformaldehyde in 0.1 M phosphate buffer saline (PBS) was perfused via the heart to fix the brain. The extracted brain tissues were post-fixed as follows: the tissues were immersed in 4% paraformaldehyde in 0.1 M PBS for 24 hours, then immersed in 10% paraformaldehyde in 0.1 M PBS for 1 week, and then sequentially immersed in 5%, 10%, and 20% sucrose in 0.1 M PBS for 1 week each.
The region containing the IONs was sectioned from the immersed brain tissue and 30-µm-thick coronal slices were prepared. The brain sections were stained with Nissl stain.

Experimental environment
The experimental environment was constructed in our previous research 1 . To encourage the rats to keep standing, the water supply was placed at a height that required the rats to stand bipedally, and water was continuously provided to reward rats for standing. The force between the water supply and mouth of the rat was monitored with a three-dimensional force sensor (Tec Gihan, Japan). The tail of each rat was tethered to the environment with a string and connected with a load cell (Minebea, Japan) to prevent its use to support the bipedal posture and to monitor the force generated at the tail. A previous study of intact rats showed that, at maximum, about 50% and 10% of the torque required for stabilization was generated from the mouth and tail, respectively 1 . A force plate (TF-2020-A; Tec Gihan) was placed on the floor of the standing rat and the floor reaction force was measured at a measurement frequency of 1000 Hz.
Measurement of the motion A motion capture system (Oqus300+; Qualisys, Sweden) was placed around the rat to measure its movement. Reflective markers were attached to the rat's skin overlying 21 body landmarks. The markers at the center of the body were as follows: the top of the head (H1), the centers of the right and left scapulae (ShM), one-and two-thirds of the H1 and ShM (H2 and H3), and the center of the right and left iliac crests (HipM). Furthermore, the following markers were placed on the right and left sides of the body: scapulae (ShR and ShL), elbow joint (ElR and ElL), distal head of the ulna (ArmR and ArmL), iliac crest (HipR and HipL), greater trochanter (GtR and GtL), lateral condyle of the knee (KneeR and KneeL), lateral malleolus (AnkleR and AnkleL), and the distal end of the fifth metatarsal (MtR and MtL).
The measurement frequency of the motion capture system was 500 Hz.

Analysis of intersegmental coordination Intersegmental coordination was analyzed using similar methods
to those of previous studies on standing intact rats 1 and standing humans 3 . The time series of the elevation angles of the four body segments (trunk, thigh, shank, and foot) on the sagittal plane were calculated from the measured positions of the head, greater trochanter, knee, heel, and metatarsus. A matrix of the segmental motion was composed by arranging the time series of the elevation angles, and singular value decomposition (SVD) was applied to this matrix to derive the principle components of the segmental motion 1 . Through this decomposition, we can obtain the combination of body segments that are simultaneously active. This is called intersegmental coordination. To investigate the effect of an ION lesion on intersegmental coordination, the component of the intersegmental coordination and the contribution rate of each intersegmental coordination to the movement were compared between intact and IO rats.

Analysis of the center of mass motion
In our previous study, we measured the ratio between the COM of each body segment to the endpoints of the segment and the weight ratio of each segment to body weight in rats 1 . These data were used to calculate the position of the COM during standing. The endpoints of the trunk, thigh, shank, and foot were derived from the motion capture data. Here, the positions of the left and right body markers were averaged. Then, using the relationship between the endpoint position and the COM, the position of the COM of each body segment was obtained. The COM of the body was calculated by the weighted sum of the COM in each body segment.
The probability density function (PDF) of the COM position in the sagittal plane was derived from the measured COM motion and compared between intact and IO rats. The PDF of COM assesses the features of the distribution of the COM similarly to the histogram of the COM. In contrast to the discrete bin-by-bin existence probability obtained by histogram, the PDF derived using kernel distribution can represent a continuous distribution of the COM. We used the "ksdensity" function in MATLAB to compute the PDF.
To show the frequency characteristics of body sway, we derived the power spectrum density (PSD) of the COM motion. The PSD was calculated using the maximum entropy method (MEM), which is stable in the estimation in the low frequency range (see Funato et al. 2017 1 for a comparison of MEM and fast Fourier transform for PSD estimation in rats). To calculate the MEM, the Burg method with 512 dimensions was used 4,5 . Because the COM movement in rats during bipedal standing had a characteristic frequency around 1 Hz 1 , the PSD was analyzed in the range of 0.05-10 Hz.

Mathematical model of posture control in rats
We model the body motion of the standing rat using an inverted pendulum with one link from the ankle to the COM as follows.
where θ is the elevation angle of the body from the vertical, m, h, J, and g are the mass, length, and moment of inertia of the pendulum and acceleration due to gravity, respectively, and is the control torque to maintain posture. We assume the following feedback control model with up to third-order nonlinear control terms (Fig. 4).
where θΔ is θ containing the sensory delay, kP, kD, and kI are the proportional, differential, and integral control gain, respectively, and σξ is the noise. The second-order nonlinear control term kP1 acts on the skewness of the COM distribution (the difference between forward and backward motion). Here, the skewness in the rat COM distribution seemed minor (Fig. 3A) and, to keep the model simple, we set kP1 = 0.  � solution: restraint condition:

Analysis of the mathematical model
Here, for the purpose of simplicity, the integral term is approximated by Substituting Eqs. (S7) and (S8) into Eq. (S6), the equation becomes and by multiplying e -jωt by both members, the above equation becomes Here, from the following relationship Eq. (S10) becomes By integrating Eq. (S12) for one cycle, the equation becomes Using the relationship − Δ = cos Δ − sin Δ , Eq. (S13) becomes We can extract a cyclic solution from this equation. From the cyclic solution a = re jφ and Eq. (S14), By comparing the real and imaginary parts of Eq. (S15), the following relationship is obtained.
Here, we consider that Δω is small, namely, sinΔω ≅ Δω and cosΔω ≅ 1, and the equation becomes From the above equations with ̇= 0, ̇= 0, the following cyclic solution is obtained.
The equation of angular frequency ω shows that the frequency of the periodic solution is determined by the differential control gain kD and the integral control gain kI. Furthermore, the integral control gain kI has a relatively smaller effect than kD because the square of the Δ ≅ 0.04 s 6 works with kI. Therefore, approximating kI ≅ 0, we obtain . (S18) The same approximation of kI ≅ 0 for the amplitude r yields (S19)

Identification of control parameters
To investigate the behavior of the mathematical model, we performed a computer simulation. A numerical calculation of the equation of motion including noise (stochastic differential equation) was performed by using the Euler-Maruyama method. The sampling frequency was was obtained by the simulation and it was converted to the position of the COM using the trigonometric function h sinθ. Then, the PDF and PSD of the simulation were calculated from the COM motion obtained, applying the same method used for the rat experiments.
The quantitative evaluation of the control system, that is, identification, was performed by a comparison of the behavior in the simulation and rat experiments. The control gains kP0, kP2, kD, and kI and noise σ of the simulated model were taken as unknown parameters, and we repeatedly performed the simulation while changing these values. Here, the delay time Δ was set as 0.04 s based on previous research 6 .
The values for the unknown parameters were determined so that the PSD and PDF of the simulation and experiment were as close as possible. The evaluation function was set as the sum of the root mean square (RMS) of the difference in the PDF and RMS of the logarithmic difference in the PSD. We searched for parameters that minimized the evaluation function. The Generic Algorithm was used for the parameter search. We used the "ga" function in MATLAB. To assess the convergence of the parameters, the identification was repeated ten times with different random seeds, and ten types of parameters were obtained for each experimental trial. Next, we tested the convergence of the identification by comparing the variation in the ten parameter sets obtained with the variation among experimental trials. Then, the optimal parameter set with the smallest evaluation function among the ten identification trials was adopted as the parameters of the model. The optimal parameter sets for intact and IO rats were compared to consider the effect of the ION lesion.