A novel IMU-based clinical assessment protocol for Axial Spondyloarthritis: a protocol validation study

Experimental data collection

Eleven participants (age 27.3 ± 5.5 years, mass 76.1 ± 7.7 kg, height 1.81 ± 0.06 m, one of which female) gave written informed consent to take part in this study and ethics were approved by the Research Ethics Approval Committee for Health (EP 17/18 128) at the University of Bath. Participants were required to be between 18 and 40 years old and have no history of bone, joint, or neuromuscular problems, or a current musculoskeletal injury.

Inertial measurement units (Avanti; Delsys inc., Natick, MA, USA) used in this study were equipped with a triaxial accelerometer (±2 g range), a triaxial gyroscope (±2000° s⁻¹ range), and a triaxial magnetometer. Wireless IMU data was recorded using Delys software (EMGworks Acquisition 4.5.4; Delsys inc., Natick, MA, USA) which included the application of a Delsys proprietary Kalman filter to convert the multi-axial sensor data into a quaternion-based orientation. The sensor data was sampled at 74 Hz and then downsampled to 50 Hz. IMU position and orientation in 3D space was compared to gold standard motion capture using 12 infrared cameras (Oqus; Qualisys, Göteborg, Sweden) and tracked via custom 3D printed clusters attached to each IMU device (Fig. 1A). Motion capture data was collected using Qualisys Track Manager (Qualisys QTM 2019.2; Qualisys AB, Göteborg, Sweden) at a sampling frequency of 100 Hz and then downsampled to 50 Hz. Both IMU and motion capture data were filtered using a zero-phase low pass 2nd order Butterworth filter with a cut-off frequency of 5 Hz (Charry, Umer & Taylor, 2011). The two systems were triggered simultaneously via a custom triggering system sending a TTL signal from the EMG station to the motion capture system.

Sensor locations were based on Lee et al. (2011b) with minor modifications. Each participant was equipped with 5 IMUs: one sensor was attached using double-sided tape directly on the skin (Fig. 1B) at the height of the spinous processes of sacrum one (S1). Three sensors were attached to foam supports (Fig. 1B) that acted to reduce soft tissue artefact by increasing the attachment surface with the skin. These sensors were placed over first lumbar vertebra (L1), first and sixth thoracic vertebra (T1 and T6). Finally, the last sensor was fixed on the head occiput with a head band. The spinous processes were located by palpation and measuring the distance between the spinous processes C7 and S1 with a flexible ruler (Ernst et al., 2013).

The protocol included five functional movements inspired by the BASMI protocol: trunk flexion, trunk extension, trunk lateral flexion, cervical rotation and cervical flexion/extension (Fig. 2). Due to technical contingencies, cervical flexion/extension was collected on five out of eleven participants. Only lumbar, lower thoracic and upper thoracic motions (sensors S1, L1, T6, T1) were recorded for trunk movements (Figs. 2A–2C). Cervical motions were excluded from the analysis as not part of the BASMI protocol, and due to the uncontrolled compensatory head movements performed by the participants whilst performing trunk lateral bending. During the cervical movements (Figs. 2D and 2E), only the cervical spine sensors (T1 and SK) were used for the analysis (Fig. 1B).

Every movement was repeated 3 times and performed at a slow constant speed, with the goal to reach the maximum ROM. Details of the instructions given to the participants are reported in Table 1. During each trial, the first 8 s were dedicated to IMU calibration to allow the Delsys Kalman filter to stabilise the orientation from the initial state. Participants stood in an upright position with their feet shoulder-width apart and looked straight ahead during this period.

Spinal angles calculation

The following analysis includes the calculation of the trunk segment angles, which are representative of the segment orientation with respect to the global reference system, and joint angles, which express the relative orientation between two consecutive sensors (sensor orientation with respect to the next, located more caudally onto the spine). The angle representations chosen are clinical angles from Crawford, Yamaguchi & Dickman (1999), where flexion/extension angles F were used for movements in the sagittal plane, lateral flexion angles L were used for the lateral flexion, and axial rotation T for the cervical rotation.

IMU segment angles calculation

The segment angles were defined as the orientation ${}_G^{B}q$ of a spinal body segment ( $B$) with respect to the global reference system $G$ (Fig. 3B). To do this, three transformations were applied (Eq. (1)): (i) the rotation of the fixed ( $F$) reference system with respect to the global ( $G$) reference system ( ${}_G^{F}q$), (ii) the sensor ( $S$) rotation with respect to the fixed ( $F$) reference system ( ${}_F^{S}q$), and (iii) the spinal segment ( $B$) rotation with respect to the sensor ( $S$) reference system ( ${}_S^{B}q$).

The transformation ${}_S^{B}q$ is representative of the soft tissue artefact’s contribution to the measure, which was not experimentally measured in this study. For this reason, the sensor ( $S$) and body ( $B$) were hypothesised being a rigid body with a constant relative orientation ( ${}_S^{B}q$), and aligned reference systems (Fig. 3B). The term ${}_S^{B}q$ is therefore a quaternion with no rotation (Eq. (2)), and the final ${}_G^{B}q$ calculation can be simplified as follows (Eq. (3)): (2) $${}_S^{B}q=\left[1,0,0,0\right]$$ (3) $${}_G^{B}q={}_G^{S}q={}_G^{F}q{}_F^{S}q$$

The sensor orientation with respect to its fixed reference system ( ${}_F^{S}q$), was calculated internally by the IMU software. The magnetometer provided the local north direction for the heading calculation (yaw), the accelerometer gravity vector contributed to the tilt estimation (roll and pitch) and the gyroscope angular velocity completed the data for the orientation calculation through sensor fusion. The sensor reference system was set up with the $\hat{Z}$ axis parallel and opposite to the gravity vector, and the $\hat{X}$ axis aligned with the local magnetic north (Fig. 3B).

IMU kinematic constraints and alignments

In the proposed method, different kinematic constraints were applied to correct for the vertical axis angular drift and improve the estimation of trunk segment angles. The kinematic constraints were applied from the end of the eight second calibration phase ( $t_{zero}$) throughout the whole recording for the sagittal and frontal plane functional movements. More specifically, the ${}_G^{F}q$ rotation constrained the movement to lie on the global sagittal or frontal plane according to the nature of the movement (see “Appendix”). This constraint assumes that the sensors reflect the actual tilt of the spine in the sagittal and frontal planes, and that the sensor heading is always aligned with the global reference system $\hat{X}$ axis. For the movement in the transversal plane, an offset calculation was performed at $t_{zero}$, and this constant offset was subtracted throughout the whole recording (see “Appendix”). This alignment assumes that the sensors reflect the actual tilt of the spine in the sagittal and frontal planes, and that the spine axial rotation is zero at $t_{zero}$.

An optimisation analysis was also conducted to further improve the estimation of the spinal segment and joint angles during lateral flexion. For each sensor, four different auxiliary axes (Fig. 3A), in addition to the $\hat{z}$ axis, were created on the sensor sagittal plane ( $\hat{x}\hat{z}$ plane). The optimisation identified the axis that generated the constraint producing the smallest measurement error, and the same axis was used for all participants (see “Appendix”).

Motion capture segment angles calibration

A set of marker clusters were attached to the IMUs to compare the alignment between their local reference system with the IMU sensors local reference system (Fig. 1A). A virtual rigid body was created for each cluster in the Qualisys Track Manager (QTM), and its orientation was exported in rotation matrix format to be subsequently converted into quaternions.

To compare the clusters and the IMUs segment angles, any misalignments of the participant’s sagittal plane with respect to the $\hat{X}\hat{Z}$ plane of the motion capture global reference system (Fig. 3A) was compensated. A vertical rotation was imposed at the instant $t_{zero}$ on the motion capture orientation to align the sensor on the global sagittal ( $\hat{X}\hat{Z}$) plane. Optical motion capture systems cannot drift and therefore kinematic constraints were not necessary. The same process was used for all the functional movements, including the calculation of the spine segment angle.

Motion capture and IMU joint angles

The joint angles of the spinal segments were obtained through a quaternion multiplication (Eq. 4).

(4) $${}_S^S{q}_{2-1}={}_G^S{\overline{q}}_1{}_G^S{q}_2$$where ${}_G^S{\overline{q}}_1$ indicates the conjugate of the unit quaternion expressing the rotation of sensor 1 with respect to the global reference system, ${}_G^S{q}_2$ indicates the unit quaternion expressing the rotation of sensor 2 with respect to the global reference system, and ${}_S^S{q}_{2-1}$ is the quaternion expressing the relative orientation of sensor 2 with respect to sensor 1 (Fig. 3B).

Flexion/extension, lateral flexion and axial rotation for segment and joint angles

Finally, by using the convention from Crawford, Yamaguchi & Dickman (1999), the segment and joint angles were expressed in quaternions and transformed into clinical angles.

(5) $$\theta ={\tan }^{-1}\left({\displaystyle \frac{{\hat{x}}_y}{{\hat{x}}_x}}\right)\varphi ={\tan }^{-1}\left({\displaystyle \frac{{\hat{x}}_y\sin \theta +{\hat{x}}_x\cos \theta }{{\hat{x}}_z}}\right)$$ (6) $$F=\varphi \cos \theta L=-\varphi \sin \theta$$ (7) $$T=-{\tan }^{-1}\left({\displaystyle \frac{-{\hat{z}}_z\sin \theta -{\hat{y}}_z\cos \theta }{{\hat{z}}_z\cos \theta -{\hat{y}}_z\sin \theta }}\right)$$where F is flexion/extension angle (Eq. (6), left), L is lateral flexion angle (Eq. (6), right), and T is axial rotation angle (Eq. (7)). An example of such angle traces is shown in Fig. 4A. It should be noted that ${\hat{z}}_z$ has an opposite sign with respect to the original formulation because that axis in the reference system of this study is directed backwords.

Range of motion and continuous spinal kinematics data segmentation

Using a custom Matlab function (MATLAB R2017b; MathWorks inc., Natick, MA, USA), segment and joint angles were processed, and spinal ROM and kinematics were calculated. The ROM data segmentation consisted of finding the angles between the anatomical reference position and the maximum ROM, for both sensors (segment angles), and spinal segments (joint angles). The only exception was cervical flexion/extension where ROM was calculated as the angle between maximum cervical flexion and maximum cervical extension. Since each movement was repeated three times and each repetition consisted of two phases (performing the movement and then returning to the reference position), flexion and extension produced 6 segments (Fig. 4B), lateral flexion and cervical rotation produced 12 segments (6 left and 6 right, Fig. 4C), and cervical flexion/extension only 3 segments (Fig. 4D). The segmentation points were calculated on IMU data, and the same points were adopted for motion capture data segmentation. Joint and segment angle ROM for both IMU and motion capture were obtained by subtracting the angular value at the beginning from the angular value at the end of every segment. This produced a ROM value for every segment. Finally, the absolute value was calculated for each of these ROMs. For the kinematic data analysis, the angular traces were segmented between the onset of the first movement and the end of the last movement (see the red lines in Figs. 4B–4D) to evaluate the sensors performance across a complete functional movement.

Sample size calculation

The Bland–Altman sample size calculation was performed on data collected on six participants. The statistical power analysis for Bland–Altman (SciStat.com, MedCalc Software Ltd, Ostend, Belgium); was set considering a type I error (alpha) of 0.05 and a type II error (beta) of 0.20. The maximum difference between IMU and motion capture measurement was set to 2° for the sagittal plane, 8° for the frontal plane and 4° for the transverse plane, considering previous studies carried out with similar technology on axSpA patients (Aranda-Valera et al., 2020). The sample size calculation yielded three different sample sizes for the three anatomical planes: 195 samples for the sagittal plane, 114 sample for the frontal plane and 73 sample for the transverse plane. Each range of motion measurement was considered as an individual sample due to the intra- and inter-variability of the measure. For the same reason a Bland–Altman repeated measure correction was not applied.

Statistical analysis

The validity of spine ROM values was assessed using correlation and Bland–Altman analyses. These tests were performed to examine the correlation, median bias and limits of agreement (95% confidence interval, 1,45 IQR (interquartile range)) between the ROM calculated via the IMU and motion capture system.

The relative reliability of the spine ROM values was assessed using the intraclass correlation coefficient ICC_2,1 (Two-way random effects, absolute agreement, and single rater). ICC values between 0.6 and 0.8 represented a good level of agreement, and >0.8 represented excellent agreement.

The Kolmogorov–Smirnov normality test was conducted to assess the distribution of the segment and joint range of motion values. After the normality test, it was chosen to perform a nonparametric analysis as the data were not normally distributed.

The absolute reliability of the spine ROM values was assessed using standardised error measurement (SEM), calculated as $\mathrm{SEM}=s.d{.}_{\mathrm{pooled}}\cdot \sqrt{1-\mathrm{ICC}}$and minimum detectable change (MDC), calculated as $\mathrm{MDC}=\mathrm{SEM}\cdot 1.96\cdot \sqrt{2}$

ROM mean absolute percentage error (MAPE) was also presented in the results. For the continuous spinal kinematic analysis, the root mean square error (RMSE), maximum absolute error (MAE), and Spearman correlation coefficient (ρ) were calculated for both segment and joint angles. Finally, the distribution of absolute and percent error frequency for segment and joint angles were calculated and plotted for visual analysis.

Analysis of spinal range of motion

The spinal ROM mean and SD for all participants are reported in Tables 2 and 3 for all the trials and all the sensors. Every table contains data of both motion capture and IMU.

Table 3:

Joint angle ROM mean and standard deviation.

ROM data (mean ± SD) for motion capture (mocap) and inertial sensors (IMU) for 11 participants (5 for cervical flexion/extension). Joint angles are reported for lumbar, lower thoracic, upper thoracic, cervical, and an overall mean. The movements from left to right were flexion, extension, lateral flexion, cervical rotation, and cervical flexion/extension.

		Flex	Ext	Lat Flex	Cer Rot	Cer F/E
	mocap	50.5 ± 10.3°	18.8 ± 6.8°	26.1 ± 4.2°	–	–
	IMU	51.7 ± 10.0°	19.1 ± 7.1°	25.8 ± 4.4°	–	–
	mocap	15.6 ± 6.9°	9.0 ± 4.8°	14.1 ± 5.6°	–	–
	IMU	16.1 ± 6.9°	9.0 ± 4.7°	14.4 ± 5.3°	–	–
	mocap	20.6 ± 8.7°	11.9 ± 6.0°	5.3 ± 4.4°	–	–
	IMU	20.6 ± 8.6°	12.4 ± 5.8°	6.0 ± 4.0°	–	–
	mocap	–	–	–	67.8 ± 9.3°	93.5 ± 11.3°
	IMU	–	–	–	68.1 ± 9.0°	93.7 ± 10.7°

DOI: 10.7717/peerj.10623/table-3

The overall mean ROM measurement of segment and joint angles across all movements had a small and comparable bias (≤0.5°), whilst the overall mean limits of agreement proved to be on average 20% higher for the joint compared to segment angle ROM (Tables 4 and 5). The lumbar and thoracic ROM measured in the sagittal plane had the smallest bias (≤0.6°) and limits of agreement (≤1.9°), whilst the measurements during lateral flexion showed smaller bias (≤1.0°) and much larger limits of agreement (≤6.4°) (Table 5). The cervical rotation joint angle (transversal plane) had bias close to zero 0.0° but the limits of agreement were much larger (≤2.3°), resulting in greater variability of the data (Table 5).

Table 5:

Joint angle analysis for ROM and kinematics.

ROM data: nonparametric Bland–Altman median bias and limits of agreement (1.45 IQR). Kinematics data: mean and standard deviation of root mean square error (RMSE), maximum absolute error (MAE) and Spearman correlation coefficient (ρ). Joint angles are reported for lumbar, lower thoracic, upper thoracic, cervical, and an overall mean. Reliability Analysis: intraclass correlation coefficient (ICC), standardised error measurement (SEM), and minimum detectable change (MDC). The movements from left to right were flexion, extension, lateral flexion, cervical rotation, and cervical flexion/extension.

		Flex	Ext	Lat Flex	Cer Rot	Cer F/E
	bias ± LOA	1.0 ± 1.7°	0.3 ± 1.6°	−1.0 ± 6.4°	–	–
	RMSE	1.9 ± 0.4°	1.4 ± 0.7°	2.4 ± 1.0°	–	–
	MAE	4.4 ± 1.0°	2.9 ± 1.6°	5.4 ± 1.8°	–	–
	ρ	0.99 ± 0.01	0.99 ± 0.01	0.99 ± 0.01	–	–
	bias ± LOA	0.3 ± 1.5°	0.1 ± 1.2°	0.3 ± 3.6°	–	–
	RMSE	1.3 ± 0.6°	1.4 ± 0.5°	1.5 ± 0.5°	–	–
	MAE	2.6 ± 0.8°	2.8 ± 1.1°	3.8 ± 1.2°	–	–
	ρ	0.98 ± 0.03	0.97 ± 0.04	0.98 ± 0.03	–	–
	bias ± LOA	0.0 ± 1.9°	0.6 ± 1.3°	0.8 ± 5.6°	–	–
	RMSE	2.7 ± 1.2°	2.8 ± 1.3°	2.9 ± 1.3°	–	–
	MAE	4.6 ± 1.3°	4.2 ± 1.6°	6.5 ± 2.2°	–	–
	ρ	0.98 ± 0.01	0.98 ± 0.03	0.48 ± 0.53	–	–
	bias ± LOA	–	–	–	0.0 ± 2.3°	0.4 ± 1.9°
	RMSE	–	–	–	1.7 ± 0.8°	2.6 ± 0.9°
	MAE	–	–	–	3.7 ± 1.2°	6.1 ± 1.0°
	ρ	–	–	–	0.92 ± 0.09	0.99 ± 0.01
Mean	bias ± LOA	0.5 ± 1.7°	0.3 ± 1.5°	0.2 ± 4.8°	0.0 ± 2.3°	0.4 ± 1.9°
	RMSE	2.0 ± 1.0°	1.8 ± 1.1°	2.3 ± 1.1°	1.7 ± 0.8°	2.6 ± 0.9°
	MAE	3.9 ± 1.3°	3.3 ± 1.6°	5.2 ± 2.0°	3.7 ± 1.2°	6.1 ± 1.0°
	ρ	0.99 ± 0.02	0.98 ± 0.03	0.82 ± 0.38	0.92 ± 0.09	0.99 ± 0.01
RA	ICC	0.998	0.989	0.952	0.986	0.993
	ICC(LB UB)	(0.997 0.998)	(0.985 0.992)	(0.942 0.961)	(0.981 0.990)	(0.981 0.998)
	SEM	0.8°	0.7°	2.1°	1.1°	0.9°
	MDC	2.1°	2.0°	5.7°	2.9°	2.4°

DOI: 10.7717/peerj.10623/table-5

The calculation of spinal ROM in the sagittal plane showed the lowest measurement error, followed by the cervical rotation and lateral flexion ROMs (Tables 4 and 5). The highest overall measurement error was observed in the frontal plane (lateral flexion) for both segment (0.5 ± 2.4°) and joint (0.2 ± 4.8°) angles, and in cervical flexion/extension segment angles (0.5 ± 2.5°) (Tables 4 and 5). This trend was confirmed in the Bland–Altman analysis, where 100% of the segment flexion (0.8% MAPE), 97% of segment extension (2.6% MAPE), and 96% of segment lateral flexion (3.5% MAPE) fell within the 10% threshold. Comparable data distribution (within 10% threshold) and MAPE values were found for the joint angles, but with higher percent errors, especially for the lateral flexion, which led to fewer data (40%) lying within the 10% percent error threshold (Fig. 5C). In fact, during the lateral flexion, the limits of agreement were relatively high across all spinal joint angles, reaching a maximum of 6.4° for the lumbar segment (Table 5), whereas the lateral flexion segment angles showed a maximum limit of agreement of 2.3° (Table 4).

The relative reliability analysis (ICC_2,1) showed excellent agreement for both segmental (Table 4) and joint angles (Table 5). The absolute reliability analysis resulted in MDC values ranging between 1.3° and 2.9° for the segmental angles, and 2.1° and 5.7° for joint angles. In both cases, spinal flexion and extension movement showed the lowest MDC values whilst lateral flexion had the highest MDC values.

Finally, the correlation analysis on the ROM data showed that the IMUs measurements were highly and significantly correlated with motion capture measurements (on average r > 0.97 for both segment and joint angles), with the lowest being the lateral flexion joint angle (r = 0.95, p < 0.001).

A Kolmogorov–Smirnov normality test was conducted on both segment and joint angles for all sagittal, frontal and transversal planes (Table 6). The transversal plane did not pass the normality test.

Kinematics analysis

Segment angle absolute error was less than 4° for more than 90% of the traces across all the movements, with the exception of cervical rotation, for which was about 85%. More specifically, the sagittal plane angle measures were the most accurate, showing an absolute error smaller than 2° and percent error smaller than 10% for more than 66% of the data (Figs. 6A and 6B). Similarly, all measurements of joint angle traces, excluding cervical rotation, had an absolute error lower than 4° for 87% of the data (Fig. 6C). The lateral flexion analysis presented two different pictures for absolute errors and percent errors. Absolute error frequency showed a constant decrease in classes containing greater errors for both segments and joint angles (Figs. 6A and 6C), whereas the percent error frequency was more evenly distributed, with lateral flexion and cervical rotation even showing an increase for the class err > 40% (Figs. 6B and 4D). On average, with cervical rotation excluded, the correlation coefficient between IMUs and motion capture segment angle traces (Table 4) was very high (ρ = 0.99). The RMSE was low for all the sensors, for all the movements (~1.5°), with 2.5° being the highest value recorded for segment T6 during flexion and extension (Table 4). There was high variability between sensors and movements, with an MAE ranging from 1.1° to 5.4°, also MAE values were higher for movements with greater ROM (Table 4).

The joint angles traces (Table 5) measured via IMUs had very high correlation coefficients compared with the motion capture analogous, the only exception being upper thoracic for lateral flexion (ρ = 0.48). The RMSE was consistently low across all spinal segments with the highest value recorded on upper thoracic during lateral flexion (2.9 ± 1.3°). Similarly, the MAE for all joint angles traces was low, except for the cervical flexion/extension angle (6.1 ± 1.0°) and the lateral flexion angle which was consistently high across all spinal segments (5.2 ± 2.0°).

Flexion-extension constraint (trunk flexion, trunk extension and cervical flexion/extension)

The $\hat{y}$ axis unit vector projection on the $\hat{X}\hat{Y}$ plane was found, and its orientation $\alpha$ with respect to the $\hat{X}$ axis was calculated through the four-quadrant inverse tangent (Eq. (A1)). This was done for all the orientations at all the time points $t$. Successively, the angles ${\beta }_t$ between the projection and the $\hat{Y}$ axis were calculated by subtracting $\pi /2$ from ${\alpha }_t$ (Eq. (A2)).

(A1) $${\alpha }_t={\tan }^{-1}\left({\hat{y}}_y\left(t\right)/{\hat{y}}_x\left(t\right)\right)$$ (A2) $${\beta }_t={\alpha }_t-\pi /2$$

Then, the quaternion angle by which the sensor orientation needs to be rotated around the fixed axis $\hat{Z}$ was: (A3) $${}_G^{F}q=q_{\beta Zt}=\left[\cos \left(-{\displaystyle \frac{{\beta }_t}{2}}\right),0,0,\sin \left(-{\displaystyle \frac{{\beta }_t}{2}}\right)\right]$$

Finally, the reorientation was performed by multiplying every quaternion ${}_G^{F}q$ by the corresponding sensor orientation quaternion ${}_F^{S}q$ (Eq. (A4)).

Lateral flexion constraints (trunk lateral flexion)

The angle $\gamma$ between the local $\hat{z}$ axis vector and its projection on the $\hat{X}\hat{Y}$ plane was calculated at $t_{zero}$ (Eq. (A5)), then a unit vector ${\hat{h}}_0=\left[001\right]$ was rotated around its $\hat{y}$ axis by an angle of $-\gamma$ (Eq. (A6)).

The quaternion rotation ${}_F^{S}q$, representing the sensor orientation in space with respect to its fixed reference system, was then applied to ${\hat{h}}_{0,-\gamma }$. The resulting unit triad has at $t_{zero}=8$ s its $\hat{z}$ axis aligned with the $\hat{X}\hat{Y}$ plane, that is now called ${\hat{h}}_t$ (Eq. (A7)).

The orientation ${\alpha }_t$ of the projection of vector $\hat{h}\left(t\right)$ on the $\hat{X}\hat{Y}$ plane was then calculated (Eq. (A8)). This was done for all orientations at all time points $t$. Successively, the angles ${\beta }_t$ between the projection and the negative $\hat{X}$ axis were calculated by subtracting $\pi$ to ${\alpha }_t$ (Eq. (A9)).

Then, the quaternion angle by which the sensor orientation needs to be rotated around the fixed axis $\hat{Z}$ was (Eq. (A10)): (A10) $${}_G^{F}q=q_{\beta Zt}=\left[\cos \left(-{\displaystyle \frac{{\beta }_t}{2}}\right),0,0,\sin \left(-{\displaystyle \frac{{\beta }_t}{2}}\right)\right]$$

Finally, the reorientation was performed by multiplying every quaternion ${}_G^{F}q$ by the corresponding sensor orientation quaternion ${}_F^{S}q$ (Eq. (A16)).

The procedure above describes the case in which the horizontal axis ( $\hat{h}$ axis) is adopted as a constraint. To adopt any of the other 4 axes ( ${\hat{z}}_2$, $\hat{z}$, ${\hat{z}}_{-1}$, ${\hat{z}}_{-2}$) as a constraint, Eqs. (A6) and (A7) become respectively: Eq. (A12) left and right (case ${\hat{z}}_2$ axis), Eq. (A13) left and right (case $\hat{z}$ axis), Eq. (A14) left and right (case ${\hat{z}}_{-1}$ axis), Eq. (A15) left and right (case ${\hat{z}}_{-2}$ axis).

Constraint optimisation (trunk lateral flexion)

For the kinematic constraints application in the lateral flexion movement, the five axes $\hat{h}$, ${\hat{z}}_2$, $\hat{z}$, ${\hat{z}}_{-1}$, ${\hat{z}}_{-2}$ (Eqs. (A7), (A12)–(A15) right) were adopted as constraints (Fig. 3A), and an optimisation that minimised the ROM measurement error absolute value was performed to find the best axes for the segment angles and the best axes pairs for the joint angles. All five axes were tested for each segment angle, and considering that the joint angles represent the relative angle between two neighbouring sensors, all twenty-five (5*5 axes) combinations were tested for the joint angles. Each of these axes was tested as a kinematic constraint, for every sensor and every participant.

Mean and standard deviation were calculated for the ROM measurement error absolute values of all sensors and participants, subsequently they were plotted against the five axes (segment angles, Fig. A1A) and the 25 axes combinations (joint angles, Fig. A1B). The segment angles constraints were numbered as follows: axis ${\hat{z}}_{-2}$ was constraint 1, ${\hat{z}}_{-1}$ was constraint 2, $\hat{h}$ was constraint 3, $\hat{z}$ was constraint 4, ${\hat{z}}_2$ was constraint 5. The 25-constraint combination numbering is shown in Table A1. The first 5 combined the lower sensor axis ${\hat{z}}_{-2}$ with the upper sensor axes going from ${\hat{z}}_{-2}$ to ${\hat{z}}_2$, the combinations from 6 to 10 combined the lower sensor axis ${\hat{z}}_{-1}$ with the upper sensor axes going from ${\hat{z}}_{-2}$ to ${\hat{z}}_2$, the same scheme repeats for all 25 constraints.

Figures A1A and A1B show mean and standard deviation distribution for ROM measurement error absolute values across respectively all axes (for segment angles) and all axes combinations (for joint angles).

Table A2 summarises the mean and standard deviation optima of the ROM measurement error absolute values. The optimum points are shown with a red diamond in Figs. A1A and A1B.

Cervical rotation alignment

The $\hat{z}$ axis unit vector projection on the $\hat{X}\hat{Y}$ plane was found at $t_{zero}$, and its orientation $\alpha \left(t_{zero}\right)$ with respect to the $\hat{X}$ axis was calculated through the four-quadrant inverse tangent (Eq. (A16)). Successively, the angle $\beta \left(t_{zero}\right)$ between the projection and the negative $\hat{X}$ axis was calculated by subtracting $\pi$ to $\alpha \left(t_{zero}\right)$ (Eq. (A17)).

Then, the quaternion angle by which the sensor orientation needs to be rotated around the fixed axis $\hat{Z}$ was (Eq. (A18)): (A18) $${}_G^{F}q\left(t_{zero}\right)=q_{\beta Z}\left(t_{zero}\right)=\left[\cos \left(-{\displaystyle \frac{\beta \left(t_{zero}\right)}{2}}\right),0,0,\sin \left(-{\displaystyle \frac{\beta \left(t_{zero}\right)}{2}}\right)\right]$$

Finally, the reorientation was performed by multiplying ${}_G^{F}q\left(t_{zero}\right)$ by every sensor orientation quaternion ${}_F^{S}q$ (Eq. (A19)).

Formulae statistical analysis

In the following formulae the subscript $A$ indicates the IMU data (Avanti), $Q$ the motion capture system data (Qualisys), $r$ the movement repeat, $p$ the participant, $s$ the sensor, $n$ the total number of movement repeats, $m$ the total number of participants, $o$ the total number of sensors, $t$ the sample at a specific time, $w$ the total number of samples per recording.