## 2 Materials and Methods
### Inverse Dynamics  
The solution over time of the kinematical problem provides all the quantities needed for the inverse dynamics. In particular, with reference to the chosen generalised coordinates \(\mathbf{q}_{i}\), the equation of motion includes the mass matrix \(\mathbf{M}\), the vector of the centrifugal and Coriolis generalised force \(\mathbf{Q}_{v}\) and the external generalised force \(\mathbf{Q}_{e}\)[23] and it is written in the form of virtual work for a virtual displacement \(\delta\mathbf{q}\) by including the work of the internal forces through the constraint Jacobian \(\mathbf{C}_{\mathbf{q}}\) and the Lagrange multipliers \(\mathbf{\lambda}\) in (16).  
\[\delta\mathbf{q}^{T}\mathbf{\left(M\ddot{\mathbf{q}}-\mathbf{Q}_{v}-\mathbf{Q}_{e}+\mathbf{C}_{\mathbf{q}} ^{T}\mathbf{\lambda}\right)}=0 \tag{16}\]  
Since the work of the internal forces is considered, the equation of motion (22) is satisfied for any virtual displacement \(\delta\mathbf{q}\): in order to refer the bodies' orientations to the local cartesian angular coordinates \(\overline{\mathbf{\varphi}}_{i}\), instead of the quaternion \(\mathbf{\theta}_{i}\), a coordinate change is performed through the matrix \(\mathbf{J}_{q}\)[23], in (17).  
\[\delta\mathbf{q}_{i}=\left[\begin{array}{c}\delta\mathbf{t}_{i}\\ \delta\mathbf{\theta}_{i}\end{array}\right]=\left[\begin{array}{cc}\mathbf{I}& \mathbf{0}\\ \mathbf{0}&\frac{1}{4}\mathbf{C}_{i}^{T}\end{array}\right]\left[\begin{array}{c} \delta\mathbf{t}_{i}\\ \delta\overline{\mathbf{\varphi}}_{i}\end{array}\right]=\mathbf{J}_{\mathbf{q},\dot{\mathbf{ \omega}}}\overline{\mathbf{q}}_{i}\ \rightarrow\ \delta\overline{\mathbf{q}}^{T}\mathbf{J}_{\mathbf{q}}^{T}\mathbf{\left(M\ddot{\mathbf{q}}-\mathbf{Q}_{ v}-\mathbf{Q}_{e}+\mathbf{C}_{\mathbf{q}}^{T}\mathbf{\lambda}\right)}=0 \tag{17}\]  
The generalised force vector \(\mathbf{Q}_{c}\) and the new constraint Jacobian \(\mathbf{C}_{\overline{q}}\) are defined in (18), leading to the final form of the multibody equation of motion.  
\[\left\{\begin{array}{cc}\mathbf{Q}_{c}=\mathbf{J}_{\mathbf{q}}^{T}\mathbf{\left(M\ddot{\mathbf{ q}}-\mathbf{Q}_{v}-\mathbf{Q}_{e}\right)}\\ \mathbf{C}_{\overline{q}}=\mathbf{C}_{\mathbf{q}}\mathbf{J}_{\mathbf{q}}\end{array}\right.\ \rightarrow\ \mathbf{C}_{\overline{q}}^{T}\mathbf{\lambda}+\mathbf{Q}_{c}=\mathbf{0} \tag{18}\]  
The coordinate transformation in (17) allows neglecting the part of the Jacobian associated with the constraints needed to keep the quaternions' unitary norm: then Equation (18) becomes a closed linear system in which the unknowns are the Lagrange multipliers associated with the rheonomic constraints \(\mathbf{\lambda}_{r}\), representing the driving forces associated to each degree of freedom, and the ones related to the scleronomic constraints \(\mathbf{\lambda}_{s}\)[25], as showed in (19).  
\[\left\{\begin{array}{cc}\mathbf{C}_{\overline{q}}=\left[\begin{array}{c}\mathbf{C }_{\overline{q},r}\\ \mathbf{C}_{\overline{q},s}\end{array}\right]\\ \mathbf{\lambda}=\left[\begin{array}{c}\mathbf{\lambda}_{r}\\ \mathbf{\lambda}_{s}\end{array}\right]\end{array}\right.\ \rightarrow\ \left[\begin{array}{cc}\mathbf{C}_{\overline{q},r}^{T}&\mathbf{C}_{\overline{q},s}^ {T}\end{array}\right]\left[\begin{array}{c}\mathbf{\lambda}_{r}\\ \mathbf{\lambda}_{s}\end{array}\right]=-\mathbf{Q}_{c} \tag{19}\]The scleronomic multipliers \(\Lambda_{s}\) are not equal to the joint reaction forces, because the actuation of the musculoskeletal system doesn't associate a single actuator to each degree of freedom: there are more muscles than degrees of freedom (redundancy) and, moreover, the muscular action is composed of a passive part which has to be considered as a known external force. In order to evaluate muscle forces and joint reactions, a muscle model is needed.