## 2 Materials and Methods
### Calf Muscle Reflex Model and Parameter Calculation  
In this paper, the following muscle reflex model is adopted:  
\[a(t)=a_{0}+p(L_{CE}(t-t_{0})-L_{CE0}(t-t_{0}))+dV_{CE}(t-t_{0}) \tag{1}\]  
where \(p\) and \(d\) are, respectively, the gains for muscle spindle length change and the rate of the change in length. \(L_{CE0}\) indicates the optimal fiber length. \(a_{0}\) denotes the initial activation level. \(L_{CE}\) is the current fiber length. \(V_{CE}\) is the current rate of change of the fiber length, obtained from the first differentiation of \(L_{CE}\)-\(L_{CE0}\). \(t_{0}\) represents the combined time delay caused by neural transmission and muscle electromechanical delays. Recent research has demonstrated that the muscle spindle primary afferents in passive muscle fire in direct relationship to muscle force-related variables [37]. We assume that there is no conflict between Equation (1) and the force sensing function. The parameters \(p\) and \(d\) can be alternatively regarded as stiffness and damping, and they connect force information with the deformation of the muscle fibers. A detailed discussion can be found in Section 4.2. Figure 1: Schematic of upright standing push-recovery. m is the mass of the body (excluding the foot part); \(\theta\) is the ankle joint angle; the clockwise direction is defined as positive; CoP is the point of center of pressure; GRFz is the z-axis ground reaction force; \(F_{m}\) is the calf muscle force. Figure 2: Ankle joint neural-muscular control model. Muscle fiber parameters \(L_{CE}\) and \(V_{CE}\) are converted from \(\theta\) and \(\dot{\theta}\) via ‘Ankle anatomy structure’. Based on the discovery noted in [38], the length of calf fibers increases linearly with the ankle joint angle in the range of the experiment. Thus, \(L_{CE}\) is simplified as:  
\[L_{CE}=c_{l}\theta+L_{CE0} \tag{2}\]  
where \(c_{l}\) represents this linear relationship. Substituting Equation (2) into Equation (1), the following equation is obtained:  
\[a(t)=a_{0}+pc_{l}\theta(t-t_{0})+dc_{l}\dot{\theta}(t-t_{0})=a_{0}+p_{a}\theta (t-t_{0})+d_{a}\dot{\theta}(t-t_{0}) \tag{3}\]  
For simplicity and without loss of generality, \(pc\) and \(dc\) are, respectively, represented by \(p_{a}\) and \(d_{a}\) in the rest of this study. It is appropriate to add an acceleration item in Equation (1), as acceleration has been proved to be useful and important in the literature. However, the results are not significantly improved in our experiment condition. A further discussion can be found in Section 4.1. Data (\(a\), \(\theta\) and \(\dot{\theta}\)) recorded from experiments are used to calculate \(p_{a}\) and \(d_{a}\). As high-level command may be involved in the control process, it is uncertain what length of time interval is most appropriate to be used. A time-increasing searching method, as illustrated in Figure 3, is proposed to find a proper time interval. The start-point of the time interval is at the time point where angular velocity begins to increase. The threshold of angular velocity is empirically set to be 0.01 rad/s. The end-point moves from the start-point to the end-point, at which muscle activation returns to the normal level with a time increment of 10 ms.  
\(p_{a}\) and \(d_{a}\) are calculated by data (\(a\), \(\theta\) and \(\dot{\theta}\)) recorded from experiments using a simple linear regress algorithm. The input is \(\mathbf{X}\) (\(\mathbf{X}\) = [\(\theta\)\(\dot{\theta}\)\(1\)]T) and the output is \(\mathbf{a}\). The parameter vector (\(\mathbf{w}\) = [\(p_{a}\)\(d_{a}\) \(a_{0}\)]) is calculated by (\(\mathbf{X}\)T)\({}^{-1}\mathbf{X}\)T\({}^{a}\). \(p_{a}\) and \(d_{a}\) are calculated at each time increment. As a result, a set of (\(p_{a}\), \(d_{a}\)) pairs is obtained.