## 2 Materials and Methods
### Static Optimization  
Once the muscle path is defined, the muscle forces have to be turned into generalised forces through the muscle Jacobian \(\mathbf{\Phi}_{q}\) by the principle of virtual work in order to include them in the equation of motion [17; 18; 25]. In general, each muscle is modelled as a series of straight-line and curved musculotendon actuators (Figure 6), so the muscle's virtual work \(\delta W_{m}\) is calculated with the sum of the individual units work considering that the muscle force \(F_{m}\) is constant along the muscle path and that there is no friction between the muscle path and the wrapping surfaces [35]. Figure 6: Muscle series. With reference to Figure 6, the three bodies \(i\), \(j\) and \(k\) are linked by the same muscle, together with an intermediate wrapping surface fixed to another body \(w\): locating the musculotonon units action lines by the unit vectors \(\vec{l}\) and the position \(\mathbf{r}\) of the attachment points, the aim is to relate the virtual work to the bodies' generalised coordinates \(\delta q\) in order to build the muscle Jacobian \(\mathbf{\Phi}_{q}\) in (27). \[ \begin{cases}\begin{array}{c}\mathbf{r}_{A}=\mathbf{t}_{i}+\mathbf{R}_{i}\overline{\mathbf{u }}_{A}\\ \mathbf{r}_{B}=\mathbf{t}_{j}+\mathbf{R}_{j}\overline{\mathbf{u}}_{B}\\ \mathbf{r}_{D}=\mathbf{t}_{w}+\mathbf{R}_{w}\overline{\mathbf{u}}_{D}\\ \mathbf{r}_{\mathbb{C}}=\mathbf{t}_{k}+\mathbf{R}_{k}\overline{\mathbf{u}}_{\mathbb{C}}\end{array} \end{cases}\begin{cases}\begin{array}{c}\mathbf{\lambda}_{AB}=\frac{\mathbf{r}_{B}- \mathbf{r}_{A}}{|\mathbf{r}_{B}-\mathbf{r}_{A}|}\\ \mathbf{\lambda}_{Bp}=\frac{\mathbf{r}_{D}-\mathbf{r}_{B}}{|\mathbf{r}_{D}-\mathbf{r}_{D}|}\\ \mathbf{\lambda}_{\mathbb{C}}=\frac{\mathbf{r}_{\mathbb{C}}-\mathbf{r}_{\mathbb{C}}}{|\mathbf{ r}_{\mathbb{C}}-\mathbf{r}_{\mathbb{C}}|}\end{array}\end{cases}\begin{cases}\mathbf{ \lambda}W_{AB}=-F_{n}\mathbf{I}_{AB}^{T}(\delta\mathbf{r}_{B}-\delta\mathbf{r}_{A})\\ \mathbf{\lambda}W_{\mathbb{B}p}=-F_{n}\mathbf{I}_{Bp}^{T}(\delta\mathbf{r}_{\mathbb{C}}- \delta\mathbf{r}_{B})\\ \mathbf{\lambda}W_{\mathbb{C}\mathbb{C}}=-F_{n}\mathbf{I}_{\mathbb{C}}^{T}\big{(}\delta \mathbf{r}_{\mathbb{C}}-\delta\mathbf{r}_{\mathbb{C}}\big{)}\end{cases} \end{cases} \tag{27}\]  
Then, once the muscle Jacobian \(\mathbf{\Phi}_{q}\) and the muscle force vector \(F_{m}\) have been assembled, the generalised muscle force \(\mathbf{Q}_{m}\) can be included as an external force in the equation of motion (16) instead of the constraint Jacobian part related to the rheonomic constraints, as written in (28). \[\mathbf{Q}_{m}^{T}\delta q=F_{m}^{T}\mathbf{\Phi}_{q}\delta q\ \to\ \delta q^{T}\Big{(}\mathbf{M} \overline{\mathbf{q}}-\mathbf{Q}_{v}-\mathbf{Q}_{e}-\mathbf{\Phi}_{q}^{T}\mathbf{F}_{m}+\mathbf{C}_{ q,v}^{T}\mathbf{\lambda}_{s}+\mathbf{C}_{q,v}^{T}\mathbf{\lambda}_{b}\Big{)}=0 \tag{28}\]  
Changing the orientation coordinates through the matrix \(\mathbf{J}_{q}\), so as to neglect the constraint Jacobian part related to the unitary norm of the quaternions, Equation (29) is obtained. \[\mathbf{\Phi}_{\overline{q}}=\mathbf{\Phi}_{q}\mathbf{J}_{q}\ \to\ \mathbf{\Phi}_{\overline{q}}^{T}\mathbf{F}_{m}-\mathbf{C}_{ \overline{q},s}^{T}\mathbf{\lambda}_{s}=\mathbf{Q}_{c} \tag{29}\]  
Each muscle force \(F_{m,i}\) depends on the muscle length \(l_{m,i}\), the muscle deformation velocity \(v_{m,i}\) and the pennation angle \(\alpha_{p,i}\) through the Hill model: these quantities can be calculated through Equation (21), evaluating the musculotonon length \(l_{mt,i}\) and deformation velocity \(v_{mt,i}\). Analysing each muscle series \(i\), in general it can be composed by \(n_{s}\) straight-line units (from \(A_{j}\) to \(B_{j}\)) and \(n_{w}\) wrapped units (from \(P_{k}\) to \(Q_{k}\)), so the musculotonon length \(l_{mt,i}\) and deformation velocity \(v_{mt,i}\) can be calculated as the sum of the individual contributions in (30). \[