A Wearable System Based on Multiple Magnetic and Inertial Measurement Units for Spine Mobility Assessment: A Reliability Study for the Evaluation of Ankylosing Spondylitis

Spinal mobility assessment is essential for the diagnostic of patients with ankylosing spondylitis. BASMI is a routine clinical evaluation of the spine; its measurements are made with goniometers and tape measures, implying systematic errors, subjectivity, and low sensitivity. Therefore, it is crucial to develop better mobility assessment methods. The design, implementation, and evaluation of a novel system for assessing the entire spine’s motion are presented. It consists of 16 magnetic and inertial measurement units (MIMUs) communicated wirelessly with a computer. The system evaluates the patient’s movements by implementing a sensor fusion of the triaxial gyroscope, accelerometer, and magnetometer signals using a Kalman filter. Fifteen healthy participants were assessed with the system through six movements involving the entire spine to calculate continuous kinematics and maximum range of motion (RoM). The intrarater reliability was computed over the observed RoM, showing excellent reliability levels (intraclass correlation >0.9) in five of the six movements. The results demonstrate the feasibility of the system for further clinical studies with patients. The system has the potential to improve the BASMI method. To the best of our knowledge, our system involves the highest number of sensors, thus providing more objective information than current similar systems.


Introduction
Ankylosing spondylitis (AS) is among the 200 specific disorders comprising rheumatic diseases, which are some of the leading causes of morbidity in the world. AS mainly affects the spine, the sacroiliac joint, and entheses [1]. AS is a chronic degenerative disease that begins with inflammatory back pain and severe inflammation of the sacroiliac joints. In more advanced stages, it can also affect the cervical spine. Chronic inflammation of the axial joints generates bony bridges, called syndesmophytes, which can join two vertebrae, limiting the movement of the spine and, in more severe states, limiting lung capacity [2,3].
The prevalence of AS in the Mexican population is around 0.6% to 0.9% [4]. Burgos-Vargas et al. [5] reported that the incidence of new cases of AS in their study of the Mexican population occurred mainly in young adults between 15 and 30 years. Despite the prevalence and the importance of the affected musculoskeletal structures, its diagnosis can be delayed by 3-11 years [6]; the Assessment of Spondyloarthritis International Society (ASAS) has designed evaluation and diagnostic criteria for patients with high suspicion [7]. The current medical tests include (a) blood tests, (b) medical imaging (CAT and MRI), (c) inflammatory activity (BASDAI index), (d) functional capacity (BASFI index), and (e) metrology (Bath Ankylosing Spondylitis Metrology Index "BASMI"). Mobility assessment of the spinal joint and the sacroiliac articulations is essential for diagnosis and evaluation of the patient's state. However, medical experts apply the BASMI index using a measuring tape and goniometers to obtain the measurements. Thus, systematic and subjective errors can appear as a function of medical experience, the correct use of the instruments mentioned above, erratic or compensatory movements of patients, or observation errors. All of these factors cause the BASMI measurements to lack accuracy and repeatability, according to some studies of the current clinical assessment [9][10][11].
In addition, the use of measuring tapes and goniometers makes the BASMI lose sensitivity to changes in the patient's status and progress; on the other hand, the spine and the sacroiliac joints generate complex movements, which cannot be evaluated with the BASMI method. Therefore, it is important to research and develop new technologybased mobility assessment methods to extensively assess the joints directly with adequate precision and repeatability, along with sensitivity to changes in information.
Optical marker-based systems have been the most widely used for tracking motion and evaluating spinal mobility [12,13]. Garrido-Castro et al. [14] established a new index (UCOASMI) correlating BASMI with optical measurements. Even though these systems have shown great accuracy and reliability, they are expensive, viable only indoors, and need a complex setup with external sources such as cameras and spaces explicitly conditioned for these purposes, limiting their use in standard clinical conditions.
MIMU-based systems for spine mobility assessment present some advantages; they are completely self-contained, allow ambulatory tracking, and can be unobtrusively attached to the patient's spine. Unlike traditional assessing methods, MIMUs can measure the continuous variation of angles in dynamic movements and maximum range of motion [27,28]. Some works studied the accuracy of MIMU-based systems, correlating them with optical systems; they showed a high level of agreement between methods and greater feasibility of MIMU systems in clinical practice [29][30][31].
Recently, some studies assessed AS patients using MIMU systems. Aranda-Varela et al. [31] established the IUCOASMI index, correlating the BASMI with their obtained metrics (flexion/extension, lateral flexion, and rotation at lumbar and cervical region) using two ViMove brand inertial sensors. O'Grady et al. [32] used two ViMove sensors to evaluate the reliability in measuring spinal range of motion. They studied the trunk flexion/extension, trunk lateral flexion, and rotation (not including measures of cervical mobility) under supervised and unsupervised conditions; they also calculated a normalized index IMU-ASMI and evaluated its reliability. Franco et al. [33] evaluated five movements, trunk flexion, trunk extension, trunk lateral flexion, cervical rotation, and cervical flexion/extension, using five commercial inertial sensors. Gardiner et al. [34] evaluated the maximum range of motion at the cervical and lumbar spine with two ViMove brand inertial sensors. In previous works [35,36], the authors presented a preliminary test with only six sensors and evaluated flexion trunk and lateral flexion in two AS patients.
Even though the studies mentioned above provide accuracy and repeatability to the measurements, they are limited by the number of sensors they use and, therefore, the amount of information they can provide for subtle and complex movements. With this motivation, we developed a novel MIMU-based system for the evaluation of spine movements. The system consists of an array of 16 small MIMU units specifically designed to be placed along the entire spine ( Figure 2a). We present the results of an experimental study aimed at validating the reliability and repeatability of the system. To the best of our knowledge, this is the only MIMU-based system capable of simultaneously monitoring different movements of the lumbar, thoracic, and cervical areas, considering the spine's serial kinematic chain with multiple degrees of freedom in various configurations. Our results indicate that the system is reliable and feasible for clinical studies and applications. It paves the way for the design of better clinical evaluation methods, with objective metrics of spine mobility, especially in rheumatology and orthopedics.

Multi-MIMU System
The proposed system's architecture consists of 16 small magnetic and inertial units (Invensense MPU-9250) and a wireless acquisition control unit. A printed circuit board was designed with the aim of having a sensor with small enough dimensions (11.7 mm × 9.3 mm) suitable for monitoring subtle movements without limiting mobility. In this way, an array of multiple sensors can be easily placed in the subjects' back and head to provide kinetic information of the entire spine at the cervical, thoracic, and lumbar sites.
Each sensor contains a triaxial accelerometer, a triaxial gyroscope, and a triaxial magnetometer (±2 g range, ±250 • /s range, ±4800 µT full range, respectively) to evaluate the subject's movements in three dimensions. The signals of the sensors were recollected through a compact and portable wireless control unit mounted in the participant's hip with a harness. The harness was designed to be comfortably worn in soft and hygienic material without limiting mobility during the assessment (Figure 2a).
The control unit consists of a digital signal processor (DSP-Teensy 3.2), interfaced with the sensors through two eight-channel fast multiplexors (TCA9548A) for I2C communication. The DSP collects the raw data from the sensors in real time and sends the information to a computer through a Bluetooth wireless communication protocol. A software desktop application was implemented to concurrently receive the raw data of the sensors and then estimate the orientation of each sensor, using a set of instances of a Kalman filter algorithm, with each instance specifically associated and calibrated per sensor [37]. The complete solution provides kinematic pose estimation of the participant's spinal mobility at a frame rate of 20 Hz per sensor ( Figure 3).

Movement Estimation Algorithm
The system comprises 16 MIMUs communicating wirelessly with a computer via a Bluetooth protocol. Each MIMU contains a triaxial accelerometer, gyroscope, and magnetometer. They are placed on different spine segments to evaluate movements in three dimensions. The tilt angles (roll and pitch) and the heading angle (yaw) are calculated through the numerical integration of the angular velocity provided by the gyroscope in the y-, z-, and x-axis according to the reference frame ( Figure 2b). It is known that the gyroscopes present a drift problem due to the accumulative errors during the recursive calculus.
The problem is solved through the implementation of a sensor fusion of the gyroscope, accelerometer, and magnetometer signals using a Kalman filter (KF) [38][39][40]. With the KF, each MIMU orientation is estimated in two phases; it is first predicted using the gyroscope signals, and then corrected by the spatial references provided by the accelerometer and The problem is solved through the implementation of a sensor fusion of gyroscope, accelerometer, and magnetometer signals using a Kalman filter (KF) [38, [40]. With the KF, each MIMU orientation is estimated in two phases; it is first predict using the gyroscope signals, and then corrected by the spatial references provided by accelerometer and magnetometer signals. The accelerometers for the tilt angles and magnetometers for the heading angle ( Figure 4).  According to the KF proposed by Lee et al. [39] and Ligorio and Sabatini [40] for estimating the two tilt angles, we implemented a modified algorithm version to include the heading angle. The purpose is to estimate the vectors S X and S Z from the rotation matrix R, which allows the coordinate transformation from the sensor frame S to the inertial frame I. On the basis of the ZYX Euler angles, R is expressed as where α is the yaw or heading angle, β is the pitch angle, and γ is the roll angle, known as tilt angles. The vector S Z (i.e., the last row of the matrix R) in Equation (1) is expressed in terms of γ and β. Hence, knowing S Z, tilt angles can be calculated as follows: Similarly, knowing S X (i.e., the first row of the matrix R) and roll and pitch angles, heading angle α can be determined by where ω is the angular velocity, and n G is the measurement noise that is assumed to be zero-mean white Gaussian. The signals from the accelerometer are modeled as where g is the gravity vector with respect to the sensor frame, defined as } × S Z where } is 9.8 m s 2 , a is the external acceleration, and n A is the measurement noise assumed as zeromean white Gaussian. As established by Luinge et al. [41], in Equation (5), the external acceleration a can be modeled as a first-order low-pass filtered white noise process.
where c a , which determines the cutoff frequency, is a dimensionless constant between 0 and 1, and ε t is the time-varying error of the acceleration model. Lastly, the signals from the magnetometers are modeled as setting h (the actual Earth's magnetic field vector) as where A −1 is the inverse of the soft-iron interference matrix, b is the hard-iron interference matrix, and h m is the distorted Earth's magnetic field (i.e., in the presence of soft-iron and hard-iron interferences) [42]. As in the gyroscope and accelerometer signals model, n M is the measurement noise assumed as zero-mean white Gaussian. It is important to develop a previous characterization step of the sensors, to set the values of the measurement noises (i.e., n G , n A , and n M ), as well as the soft-iron and hard-iron interference matrices [37,43].

Kalman Filter Design
The KF is defined by the following process model [44]: and the following measurement model: where x t in Equation (9) is the state vector, defined as A is the state transition matrix that relates the state at a previous time step t − 1 to the state at the current step t, and w is the white Gaussian process noise. In Equation (10), z t is the measurement vector, H is the observation matrix that relates the state (i.e., x t ) to the measurement z t , and v is the white Gaussian measurement noise. The minus superscript in the state vector denotes the a priori estimate. The orientation is a priori estimated by the integration of the gyroscope signals, giving the following process model: where ∆t is the time interval and ω is the skew-symmetric matrix function of the vector ω, denoted as Due to the noisy measurements of the angular velocity in Equation (4), Equation (9) must be expressed using the current gyroscope output (i.e., y G = ω + n G ) as Following the mathematical deduction reported by Lee et al. [39], Equation (13) can be developed as Thus, from Equation (14), the transition matrix A t−1 and the process noise w t−1 can be defined as Then, the process noise covariance matrix Q t−1 defined by E w t−1 w T t−1 can be redefined using Equation (16) as where Σ G , defined by E n G n T G , is the covariance matrix of the gyroscope measurements noise and is established equal to In Equation (18), σ 2 G is obtained in a previous characterization step and is a 3 × 3 matrix where the gyroscope noise variances of the x-, y-, and z-directions are in the main diagonal.
The measurement model is based on the accelerometer and magnetometer measurements since they give the spatial reference to correct the estimation error in the process model. As in [41], the error of the predicted acceleration is defined as where a − t is the a priori estimate of the external acceleration of the current time step and is defined as c a a + t−1 , which is available from a previous time. Note that the plus superscript denotes the a posteriori estimate after the filter correction.
Therefore, using Equations (6) and (19), Equation (5) can be rewritten as In Equation (18), is obtained in a previous characterization step and is a 3 × 3 matrix where the gyroscope noise variances of the x-, y-, and z-directions are in the main diagonal.
The measurement model is based on the accelerometer and magnetometer measurements since they give the spatial reference to correct the estimation error in the process model. As in [41], the error of the predicted acceleration is defined as where is the a priori estimate of the external acceleration of the current time step and is defined as , which is available from a previous time. Note that the plus superscript denotes the a posteriori estimate after the filter correction.
Therefore, using Equations (6) and (19), Equation (5) can be rewritten as From Equation (7), ℎ, the actual Earth's magnetic field vector, is where is the value of the local Earth's magnetic field. Hence, from Equations (20) and (21), the measurement model can be established sd follows: Accordingly, the measurement vector , the observation matrix , and the From Equation (7) In Equation (18), is obtained in a previous characterization step and is a 3 × 3 matrix where the gyroscope noise variances of the x-, y-, and z-directions are in the main diagonal.
The measurement model is based on the accelerometer and magnetometer measurements since they give the spatial reference to correct the estimation error in the process model. As in [41], the error of the predicted acceleration is defined as where is the a priori estimate of the external acceleration of the current time step and is defined as , which is available from a previous time. Note that the plus superscript denotes the a posteriori estimate after the filter correction.
Therefore, using Equations (6) and (19), Equation (5) can be rewritten as From Equation (7), ℎ, the actual Earth's magnetic field vector, is where is the value of the local Earth's magnetic field. Hence, from Equations (20) and (21), the measurement model can be established sd follows: Accordingly, the measurement vector , the observation matrix , and the where re is the a priori estimate of the external acceleration of the current time step and efined as , which is available from a previous time. Note that the plus erscript denotes the a posteriori estimate after the filter correction.
Therefore, using Equations (6) and (19), Equation (5) can be rewritten as From Equation (7), ℎ, the actual Earth's magnetic field vector, is re is the value of the local Earth's magnetic field. Hence, from Equations (20) and , the measurement model can be established sd follows: Accordingly, the measurement vector , the observation matrix , and the surement noise are Since , , , and are uncorrelated, the measurement noise covariance matrix, defined by , results in re , , and are the covariance matrices of the magnetometer measurement e, the acceleration model error, and the accelerometer measurement noise, ectively.
is defined as , which is the 3 × 3 matrix where the magnetometer noise iances of XYZ directions are in the main diagonal. As in [39], is defined as ‖ ‖ , where ‖ ‖ is the square of the vector norm, and is set as , the 3 × atrix where the accelerometer noise variances of XYZ directions are in the main onal. Eventually, once are estimated, the external acceleration can be ulated by Experimental Study

The Objective of the Study
The objective of the study was to assess the reliability and repeatability of the posed system for the clinical evaluation of spine motility. To this end, six movements is the value of the local Earth's magnetic field. Hence, from Equations (20) and (21), the measurement model can be established sd follows: where is the a priori estimate of the external acceleration of the current time step and is defined as , which is available from a previous time. Note that the plus superscript denotes the a posteriori estimate after the filter correction.
Therefore, using Equations (6) and (19), Equation (5) can be rewritten as From Equation (7), ℎ, the actual Earth's magnetic field vector, is ℎ = , where is the value of the local Earth's magnetic field. Hence, from Equations (20) and (21), the measurement model can be established sd follows: Accordingly, the measurement vector , the observation matrix , and the measurement noise are = , , − , = − , + .
Since , , , and are uncorrelated, the measurement noise covariance matrix, , defined by , results in where , , and are the covariance matrices of the magnetometer measurement noise, the acceleration model error, and the accelerometer measurement noise, respectively.
is defined as , which is the 3 × 3 matrix where the magnetometer noise variances of XYZ directions are in the main diagonal. As in [39], is defined as 3 ‖ ‖ , where ‖ ‖ is the square of the vector norm, and is set as , the 3 × 3 matrix where the accelerometer noise variances of XYZ directions are in the main diagonal.
Eventually, once are estimated, the external acceleration can be calculated by = , − ℊ .

The Objective of the Study
The objective of the study was to assess the reliability and repeatability of the proposed system for the clinical evaluation of spine motility. To this end, six movements The measurement model is based on the accelerometer and magnetometer measurements since they give the spatial reference to correct the estimation error in the process model. As in [41], the error of the predicted acceleration is defined as where is the a priori estimate of the external acceleration of the current time step and is defined as , which is available from a previous time. Note that the plus superscript denotes the a posteriori estimate after the filter correction.
Therefore, using Equations (6) and (19), Equation (5) can be rewritten as From Equation (7), ℎ, the actual Earth's magnetic field vector, is where is the value of the local Earth's magnetic field. Hence, from Equations (20) and (21), the measurement model can be established sd follows: Accordingly, the measurement vector , the observation matrix , and the measurement noise are = − , + .
Since , , , and are uncorrelated, the measurement noise covariance matrix, , defined by , results in where , , and are the covariance matrices of the magnetometer measurement noise, the acceleration model error, and the accelerometer measurement noise, respectively.
is defined as , which is the 3 × 3 matrix where the magnetometer noise variances of XYZ directions are in the main diagonal. As in [39], is defined as 3 ‖ ‖ , where ‖ ‖ is the square of the vector norm, and is set as , the 3 × 3 matrix where the accelerometer noise variances of XYZ directions are in the main diagonal.
Eventually, once are estimated, the external acceleration can be calculated by

The Objective of the Study
The objective of the study was to assess the reliability and repeatability of the proposed system for the clinical evaluation of spine motility. To this end, six movements Accordingly, the measurement vector z t , the observation matrix H, and the measurement noise v t are In Equation (18), is obtained in a previous characterization step and is a 3 × 3 matrix where the gyroscope noise variances of the x-, y-, and z-directions are in the main diagonal.
The measurement model is based on the accelerometer and magnetometer measurements since they give the spatial reference to correct the estimation error in the process model. As in [41], the error of the predicted acceleration is defined as where is the a priori estimate of the external acceleration of the current time step and is defined as , which is available from a previous time. Note that the plus superscript denotes the a posteriori estimate after the filter correction.
From Equation (7), ℎ, the actual Earth's magnetic field vector, is ℎ = , where is the value of the local Earth's magnetic field. Hence, from Equations (20) and (21), the measurement model can be established sd follows: Accordingly, the measurement vector , the observation matrix , and the measurement noise are = , , − , = − , + .
Since , , , and are uncorrelated, the measurement noise covariance matrix, , defined by , results in where , , and are the covariance matrices of the magnetometer measurement noise, the acceleration model error, and the accelerometer measurement noise, respectively.
is defined as , which is the 3 × 3 matrix where the magnetometer noise variances of XYZ directions are in the main diagonal. As in [39], is defined as 3 ‖ ‖ , where ‖ ‖ is the square of the vector norm, and is set as , the 3 × 3 matrix where the accelerometer noise variances of XYZ directions are in the main diagonal.
Eventually, once are estimated, the external acceleration can be calculated by = , − ℊ .

The Objective of the Study
The objective of the study was to assess the reliability and repeatability of the proposed system for the clinical evaluation of spine motility. To this end, six movements In Equation (18), is obtained in a previous characterization step and is a 3 × 3 matrix where the gyroscope noise variances of the x-, y-, and z-directions are in the main diagonal.
The measurement model is based on the accelerometer and magnetometer measurements since they give the spatial reference to correct the estimation error in the process model. As in [41], the error of the predicted acceleration is defined as where is the a priori estimate of the external acceleration of the current time step and is defined as , which is available from a previous time. Note that the plus superscript denotes the a posteriori estimate after the filter correction.
From Equation (7), ℎ, the actual Earth's magnetic field vector, is ℎ = , where is the value of the local Earth's magnetic field. Hence, from Equations (20) and (21), the measurement model can be established sd follows: Accordingly, the measurement vector , the observation matrix , and the measurement noise are = , , − , = − , + .
Since , , , and are uncorrelated, the measurement noise covariance matrix, , defined by , results in where , , and are the covariance matrices of the magnetometer measurement noise, the acceleration model error, and the accelerometer measurement noise, respectively.
is defined as , which is the 3 × 3 matrix where the magnetometer noise variances of XYZ directions are in the main diagonal. As in [39], is defined as 3 ‖ ‖ , where ‖ ‖ is the square of the vector norm, and is set as , the 3 × 3 matrix where the accelerometer noise variances of XYZ directions are in the main diagonal.
Eventually, once are estimated, the external acceleration can be calculated by = , − ℊ .

The Objective of the Study
The objective of the study was to assess the reliability and repeatability of the proposed system for the clinical evaluation of spine motility. To this end, six movements Since a − ε,t , n A , and n M are uncorrelated, the measurement noise covariance matrix, where Σ M , Σ acc , and Σ A are the covariance matrices of the magnetometer measurement noise, the acceleration model error, and the accelerometer measurement noise, respectively. Σ M is defined as σ 2 M , which is the 3 × 3 matrix where the magnetometer noise variances of XYZ directions are in the main diagonal. As in [39], Σ acc is defined as 3 −1 c 2 a a + In Equation (18), is obtained in a previous characterization step and is a 3 × 3 matrix where the gyroscope noise variances of the x-, y-, and z-directions are in the main diagonal.
The measurement model is based on the accelerometer and magnetometer measurements since they give the spatial reference to correct the estimation error in the process model. As in [41], the error of the predicted acceleration is defined as where is the a priori estimate of the external acceleration of the current time step and is defined as , which is available from a previous time. Note that the plus superscript denotes the a posteriori estimate after the filter correction.
From Equation (7), ℎ, the actual Earth's magnetic field vector, is ℎ = , where is the value of the local Earth's magnetic field. Hence, from Equations (20) and (21), the measurement model can be established sd follows: Accordingly, the measurement vector , the observation matrix , and the measurement noise are = , , − , = − , + .
Since , , , and are uncorrelated, the measurement noise covariance matrix, , defined by , results in where , , and are the covariance matrices of the magnetometer measurement noise, the acceleration model error, and the accelerometer measurement noise, respectively.
is defined as , which is the 3 × 3 matrix where the magnetometer noise

The Objective of the Study
The objective of the study was to assess the reliability and repeatability of the proposed system for the clinical evaluation of spine motility. To this end, six movements were considered: three movements involving the articulations along the thoracic and lumbar spine (anterior hip flexion, trunk lateral flexions, and trunk axial rotation) and three involving the cervical spine (cervical axial rotation, cervical flexion/extension, and cervical lateral flexion).

Recruitment
Fifteen healthy participants were invited to participate in the study (eight male and seven female). The inclusion criteria were being young adults without any mobility impairment, history of diagnosis of spine or musculoskeletal ailment, or presence of any joint or spinal pain at the time of the study. The individuals shared homogeneous anthropometric characteristics such as weight, height, and body mass index (BMI) to guarantee the same physical mobility capabilities among participants. The study was approved by the Local Ethics Committee of the General Hospital of Mexico "Dr. Eduardo Liceaga" (protocol code DI/03/17/471). All participants signed the informed consent form.

Experimental Procedure
The participants performed two repeated exercise sequences during complete sessions on 2 days each, 1 week apart from the first to the second sequence. The sessions were carried out at approximately the same period during the day on mornings. All experiments and recordings were conducted by the same operator to assess the repeatability of the system. Sociodemographic data and anthropometric measurements were collected on the first day of the study. Each day, every subject was equipped with 16 sensors along the spine, one on the head occiput to assess the cervical movements, and 15 along the back. The first back sensor was attached over the first sacral vertebra (S1) and the last was attached over the last cervical vertebra (C7), leaving an equidistant separation between each sensor. The distance between sensors was different on each subject depending on their height and natural erect posture (Figure 2a). The back sensors were mounted directly on the skin with double-sided tape, and the head sensor was mounted with a strap.
The spinal assessment protocol included six movements (Figure 2c-h), including three from the original BASMI method: cervical rotation, anterior hip flexion, and trunk side flexion. Instead of tragus-to-wall distance, we propose introducing two other head movements, the cervical flexion/extension and cervical lateral flexion. Additionally, we studied the trunk axial rotation along the lumbar and thoracic spine. The operator kindly guided the participants using standardized instructions through each movement. Participants were asked to perform their maximum effort without feeling any pain or discomfort to reach their full range of motion (ROM). Each movement was repeated three times at a slow and constant speed; before starting every series, the participants stood in a natural upright position with their feet shoulder-width apart for 5 s to set the baseline position of each sensor. All participants performed the sequence of six-movement repetitions strictly in the same order to prevent any possible bias due to arbitrary execution of the exercise routines.

Measurements
For all subjects, the Euler angles in ZYX representation were simultaneously estimated for each MIMU for all six movement series. For movement 1, movement 2, and movement 3, 15 sensor units were considered, from MIMU1 to MIMU15, placed approximately over the S1 sacral vertebra up to the C7 cervical vertebra. Thus, there were seven sensor units (MIMUs 1-7) for the lumbar spine and eight units (MIMUs 8-15) for the thoracic spine. For movement 4, movement 5, and movement 6, three sensor units were considered, from MIMU14 to MIMU16, with MIMU14 between thoracic vertebrae T1 and T2, MIMU15 at the C7 vertebra, and MIMU16 over the occipital area of the head. Then, the ranges of motion (RoM) of frontal flexions of movement 1 and cervival flexions/extensions of movement 5 were estimated from the maximum and minimum roll angles of every sensor unit. The RoM of lateral flexions of movement 2 and movement 6 were estimated from maximum and minimum pitch angles; RoM of axial rotations of movement 3 and movement 4 were estimated analogously from yaw angles.

Statistical Analysis
The uniformity of the age and anthropometric measurements was statistically evaluated using the t-test. Descriptive statistics (mean, SD, and ranges) over the RoM of the participants per MIMU, grouping for the first and second measurements, were computed.
To analyze the repeatability of the system. The intraclass correlation coefficient (ICC) was first computed for each MIMU, for each movement separately from the mean, for different pairs of measurements. Then, the grand ICC was estimated per movement, grouping the pairs of measurements of all the corresponding MIMUs involved in the respective exercise. ICC interpretation, based on [45], was as follows: <0.4 = poor, 0.41 to 0.6 = regular, 0.61 to 0.8 = good, and >0.81 = excellent. A significance level of p < 0.05 was assumed for all statistical tests. The statistical tests and analysis were carried out in Matlab R19 and SPSS 20.

Results
The study group comprised 15 healthy participants. A one-sample Student's t-test (two-tailed version) confirmed the homogeneity of participants' anthropometric characteristics, resulting in significant descriptive statistics for all the metrics, as shown in Table 1. Thirteen participants completed the entire study. Data from two participants were discarded for analysis because they had issues with the fall of some sensors Table 1. The participants' anthropometric information shows a uniform population, indicating similar physical conditions and capacities for performing the spine and head movements.

Spinal Mobility Data
Descriptive statistics of the RoM for each MIMU and every movement are reported as angles in Table 2. The RoM obtained for every section of the spine was consistent with that reported in [46][47][48]. The reliability per sensor was assessed, in general showing good to excellent correlation in movements 1, 2, 4, 5, and 6 for most MIMUs. Figure 5 presents bar plots showing the medians, ranges, and outliers obtained from each combination of 16 IMU × 2 measurements × 6 movements.
Bland-Altman plots in Figure 6 show the agreements of the RoM for every movement. As can be observed, the system presents good agreement for five of the six studied exercises, with few observed outliers, except for trunk axial rotation (Figure 6c). Table 3 presents the detailed descriptive statistics of the corresponding Bland-Altman plots in Figure 6, confirming the high repeatability for frontal and lateral flexions of the lumbar-thoracic spine and the cervical spine in axial rotation, flexion/extension, and lateral flexions. Additionally, it shows good repeatability for the axial rotation of the trunk.

Discussion
The purpose of this study was to demonstrate the reliability and repeatability of the proposed novel system. The system comprises 16 miniaturized sensors to assess the entire spine, including the cervical, thoracic, and lumbar sections, unlike previous studies [31,32,34] which used only two sensors to evaluate the cervical spine and two sensors to assess the lumbar spine in separate tests. This configuration with more sensors allows obtaining more motion information at the highest spatial resolution than existing approaches. For this reason, it is expected that the system would provide metrics more sensitive to changes in the progress of the impairment or improvements due to treatments in patients with musculoskeletal disorders, especially disabilities due to AS.
The intrarater reliability analysis for all the MIMUs (Table 2) showed mostly a good (0.61 < ICCs < 0.8) to excellent (>0.81) intraclass correlation coefficient. In particular, for the anterior hip flexion exercise (Movement 1), seven of the 15 used sensors presented a good level of reliability, as exhibited by the achieved ICC, six presented an excellent reliability level, and only two presented a moderate level. These two were placed at the bottom part of the lumbar spine.
For the trunk lateral flexion exercise (Movement 2), for the same 16 sensors, the system presented 10 sensors with a good intraclass correlation coefficient but five sensors with poor coefficients, three of them located in the lumbar region. The observed reliability level of the three sensors placed at the lumbar site showed similar behavior to the coefficient reported in [34] of the sensor located in the same anatomical region. Sensors 1 to 5 showed a poor ICC due to their position on the participants' back. For MIMU 1, a possible explanation is that it was too close to the participants' clothes, which may have induced erratic sliding of the sensor during the exercises. For sensors 2-5, the observed variability in the lumbar ICCs was likely due to anatomical factors rather than sensor error. Because participants have different body compositions (adipose tissue), a narrower channel is generated in the center of the back between the muscles, making it difficult to place the sensors and keep them in place without undesired sliding motions given the skin displacements.
For trunk axial rotation (Movement 3), the system did not give helpful information, at least not for the entire spine, in the tested configuration. The ICCs were poor in almost all sensors, except four of five sensors placed on the upper thoracic spine (MIMU [11][12][13][14][15]. A suitable explanation is that the localization of the sensors was not optimal in the selected configuration. In other words, placing them along the spine also causes erratic sliding motions of the sensors due to skin deformation during the exercises. Furthermore, a signal-to-noise problem is plausible due to the noise increment with the erratic movements of the sensors. Thus, it is necessary to continue researching the best configuration to place the sensors for axial rotation. The design presented by Molnar et al. [30] could be a better option since the sensors were placed on both sides of the spine at sites with less skin displacement. Concerning the exercises involving cervical movements (Movements 4-6), sensor 16 was attached with a strap to the head occiput. This sensor measures the range of motion of the cervical spine with respect to the back. Variability in the cervical spine ranges of motion was observed among participants, resulting in good but not excellent reliability. We suspect that not all the participants performed their maximum ranges of motion in trying to complete the series dynamically, as we were able to visually observe during the experiments. Personal natural anatomical limitations despite the anthropometric homogeneity of the population cannot be ruled out either. An option could be to record the ranges of motion in static conditions, ensuring that the maximum range has been reached, providing consistency in measurements. In any case, the dynamic recording of the exercises could give valuable temporal and differential information for future research.
As can be observed, some factors influenced the variability of the system, some attributed to the population itself, and some attributed to measurement errors. However, following the overall intraclass correlation analysis to assess the global intraobserver reliability of the system (Table 3), an excellent level of intraclass correlation was evidenced for five of the six considered movements (i.e., Movements 1, 2, 4, 5, and 6), whereas one had good ICC (Movement 3). These results indicate acceptable repeatability of the system and, thus, feasibility for further clinical studies.
From the results of the intrarater variability, it is deduced that the system offers advantages over manual clinical methods (BASMI) because it can provide the clinicians with more relevant information on how the patient performs movements, as well as how it is affected by the clinical condition, in terms of speed, smoothness of movements, continuity, and fluidity, in a more systematic and controlled fashion, thereby reducing subjectivity.
As a limitation of the current study, for the moment, we evaluated the repeatability of the system in static conditions considering the evaluation of the RoM; however, it would be interesting and require extensive research effort to assess the reliability of the system regarding the dynamic execution of the movements. This may be valuable for introducing new metrics involving temporal patterns, movements phases, and abnormal motion synergies not considered in current clinical practice.
Thus, as future work, an information reduction analysis can be done to obtain a more efficient system, with only the necessary sensors and correct locations maintained for providing important information and facilitating its practical use in clinical environments. Furthermore, the design of an ergonomic suit and belt is needed. In any case, our proposed system features the largest number of wearable sensors, allowing the investigation of new and different configurations. This can be helpful to evaluate not only the spine but also other joints such as the arms, legs, and fingers, which can also be affected by AS (or other musculoskeletal conditions). Moreover, further studies need to be conducted to confirm the inter-rater reliability involving at least two different observers.
Regarding the validity and accuracy of the system, the registered ranges of motion are consistent with the literature on spinal biomechanics [46][47][48] and similar to those reported by Franco et al. [33]. Even so, further study is needed to assess the accuracy of the system. Additionally, the study was conducted in a controlled environment where electromagnetic noise was characterized. However, for its implementation in clinical practice, it is necessary to enhance the Kalman filter to neutralize uncontrolled electromagnetic noise that may not be adequately characterized, similarly to the measurement noise covariance reported in [49,50]. Once again, more studies testing the repeatability of the instruments should be considered in the future.

Conclusions
This paper presented the development and evaluation of a novel magnetic and inertial measurement unit-based system to assess the entire spine's motion. The system comprises 16 MIMUs unlike previous studies that used at most five sensors. The applied protocol demonstrated the reliability of using the system to evaluate the entire spine, with an excellent level of intraclass correlation (>0.81) for five of the six considered movements. The present study introduced three movements (i.e., trunk rotation, cervical flexion/extension, and cervical lateral flexion), unfeasible to evaluate using manual methods such as the BASMI, the most used index for spine assessment in AS. However, it is necessary to continue researching the best configuration to place the sensors for axial rotation since the proposed approach did not give helpful information. According to the tests performed, spine evaluation using the proposed approach presents advantages over traditional methods. In addition, the quantity and configuration of the MIMUs allow obtaining more objective information on the patient's condition for a better and more timely diagnosis. The system represents a valuable tool to investigate better and more accurate evaluation metrics than current clinical methods.