Pressure based MRI-compatible muscle fascicle length and joint angle estimation

Background Functional magnetic resonance imaging (fMRI) provides critical information about the neurophysiology of the central nervous systems (CNS), posing clinical significance for the understanding of neuropathologies and advancement of rehabilitation. Typical fMRI study designs include subjects performing designed motor tasks within specific time frames, in which fMRI data are then analyzed by assuming that observed functional brain activations correspond to the designed tasks. Therefore, developing MRI-compatible sensors that enable real-time monitoring of subjects’ task performances would allow for highly accurate fMRI studies. While several MRI-compatible sensors have been developed, none have demonstrated the ability to measure individual muscle fascicle length during fMRI, which could help uncover the complexities of the peripheral and central nervous systems. Furthermore, previous MRI-compatible sensors have been focused on biologically intact populations, limiting accessibility to populations such as those who have undergone amputation. Methods We propose a lightweight, low-cost, skin impedance-insensitive pressure-based muscular motion sensor (pMMS) that provides reliable estimates of muscle fascicle length and joint angle. The muscular motions are captured through measured pressure changes in an air pocket wrapped around the muscle of interest, corresponding to its muscular motion. The muscle fascicle length and joint angle are then estimated from the measured pressure changes based on the proposed muscle-skin-sensor interaction dynamics. Furthermore, we explore an integration method of multiple pMMS systems to expand the sensor capacity of estimating muscle fascicle length and joint angle. Ultrasound imaging paired with joint encoder measurements are utilized to assess pMMS estimation accuracy of muscle fascicle length in the tibialis anterior (TA) and ankle joint angle, respectively, of five biologically intact subjects. Results We found that a single pMMS sufficiently provides robust and accurate estimations of TA muscle fascicle length and ankle joint angle during dorsiflexion at various speeds and amplitudes. Further, differential pressure readings from two pMMSs, in which each pMMS were proximally and distally placed, were able to mitigate errors due to perturbations, expanding pMMS capacity for muscle fascicle length and ankle joint angle estimation during the full range of plantar flexion and dorsiflexion. Conclusions Our results from this study demonstrate the feasibility of the pMMS system to further be incorporated in fMRI settings for real-time monitoring of subjects’ task performances, allowing sophisticated fMRI study designs.


Background
Functional magnetic resonance imaging (fMRI) is utilized for investigating brain activity in subjects performing given tasks within specific time frames during a functional scan [1]. fMRI allows for a consecutive time-series of brain activity to be compared with specific tasks and movements, which is critical data for uncovering the connections between motor dysfunction and neural activity. To synchronize fMRI frames with motor task data, and to enable real-time monitoring of subjects' task performances, various MRI-compatible sensors have been proposed to measure grip force [2], joint position [3][4][5], and net torque [6][7][8][9]. Meanwhile, the exchange of efferent motor signals from the brain and afferent sensory signals from the muscular system, through mechanoneural transduction within muscles, is known to play a critical role in fine motor control [10,11]. Therefore, the development of a MRI-compatible sensor that provides muscle fascicle length would allow for the assessment of afferent signaling within muscles. Furthermore, for populations such as persons who have undergone amputation, assessment of the subjects' task performances needs to be based on muscle fascicle movements, independent of joint kinematics and kinetics. Thus, the implementation of a MRI-safe sensor for collecting muscle fascicle data in conjunction with brain activity may be immensely important for advancing our understanding of the complex network between the central and peripheral nervous systems, and evaluating the efficacy of rehabilitation strategies to restore neuromuscular control.
Implantable sonomicrometery crystals have been reported to accurately record real time muscle fascicle length changes. Researchers have used sonomicrometery crystals to track muscle fascicle length and dynamics through direct measurement of crystal activity using a sonomicrometry amplifier [12]. Recently, it has been discovered that magnets can be used in a similar manner. Magnetic fields are able to pass through biological tissues which allows for muscle fascicle length changes to be recorded, via distance tracking between implantable magnets, without a powered system [13]. These methods have been proven successful for highly accurate muscular tracking without being affected by outside noise or skin impedance. However, both require invasive surgery and are strictly limited to a majority of the biologically intact population. Also, magnet tracking still needs to be tested for long term viability in human applications as well as MRI-compatibility.
Sonomyography via ultrasound imaging has also been utilized to derive muscle fascicle length and pennation angle measurements. This has been successful in characterizing muscle architecture for applications including gait abnormality detection, quantifying the dynamics of individual muscles, and measuring proprioceptive intent [14,15]. An advantage of ultrasound imaging is its ability to isolate muscles of interest, down to individual fascicles, resulting in precise, easy-to-visualize measurements. However, this modality is bound to clinical and laboratory settings, and the quality of the imaging is limited to trained technicians. The reliability of ultrasound imaging may also be diminished due to its reliance on individual interpretation of a given image [16].
For the purpose of avoiding the aforementioned, several efforts have focused on cross-sectional expansion of muscle during its movements [17,18]. These muscle movement-based methodologies place an air pocket, or force sensitive registers (FSRs), on the muscle of interest, and by tracking cross-sectional expansion via air pressure or contact forces, the muscle activity was estimated. Force myography (FMG) has also been explored for detecting volumetric changes of specific muscles inside a compression garment with the intention of estimating muscle activity [19]. These systems have advantages in terms of practicality for being insensitive to skin impedance changes, and do not require invasive attachment to the skin. In the previously mentioned studies, however, the general estimation strategies of variables such as muscle fascicle length and joint angles, that have non-negligible non-linearity between the sensor outputs and target variable, were not explored. In addition, these studies have neglected muscle-skin interactions and sensor dynamics, e.g. air pocket dynamics and non-linearity of FSRs.
Although most of the sensors mentioned above have proved successful in clinical settings, they remain futile in conjunction with fMRI for safety and logistical reasons. In this paper, we propose a MRI-safe, non-invasive, cost efficient, pressure-based muscular motion sensor (pMMS) and associated methodology. The pMMS accounts for muscle-skin-sensor dynamics and addresses the ability to estimate changes in muscle fascicle length and joint angles in free-space. Similar to the methodology of the previously mentioned air pocket system, the pMMS utilizes pressure changes in the air pocket due to the cross-sectional expansion of the muscle as shown in Fig. 1a. The proposed method includes a general estimation strategy of target variables that have non-negligible non-linearity with the sensor outputs, to provide robust, real-time estimates of muscle fascicle length and joint angle. These estimations from the pMMS are evaluated by comparative experiments using ultrasound imaging and a joint encoder. Through our investigation, we demonstrate the device's capability to reject perturbations from surrounding muscles by utilizing the differential pressure of two pMMSs (differential-pMMS), and demonstrate pMMSs' muscle fascicle length and joint angle estimation accuracy during the full range of motion of the ankle joint. The pMMS may provide several clinically-significant contributions including the ability to synchronize muscle fascicle length and

Mechanism and design of the pressure-based muscular movement sensor (pMMS)
The pressure-based muscular motion sensor (pMMS) consists of a small air pocket (3×0.8×0.8 cm) and a housing cell as shown in Fig. 1b. When a muscle is shortened or lengthened, the cross-sectional area of the muscle varies due to changes in the overlapping area of myosin and actin filaments [20], which are biological structures for muscular force production, and muscle fiber pennation angle [21]. By measuring pressure changes within the pMMS wrapped around a muscle of interest, the pMMS reading correlates to the cross-sectional expansion of the muscle, which allows for estimation of the target variables, muscle fascicle length and joint angle, based on the sensor dynamics as shown in Fig. 2a. The housing cell incorporates a double layer structure to minimize pressure changes due to external perturbations and isolates the muscle of interest from surrounding muscles as shown in Fig. 1.

Modeling
The schematic diagram of pMMS dynamics is shown in Fig. 2a. In this work, an internal equilibrium point x m is assumed, which has a static non-linear relationship g(.) between its variance δx m and that of the target variable δη.
x m is utilized to reflect the η-dependent steady pressure of the air pocket P a . Overall, δη results δx m based on g(.), which propagates through the muscle-skin and deforms From the ideal gas law, P a is described as where n, R, and T indicates the moles of gas in the air pocket, the gas constant, the absolute temperature, respectively, and A, x 0 , and x a are the contact area, positions of top and bottom of the air pocket, respectively. By taking the derivative of Eq.1, we obtain Because A and n were constant by the design of pMMS and T reached equilibrium with body temperature during trials, K p was considered as constant for marginal changes of δ(x 0 − x a ). The viscoelastic properties of the air pocket are known to be well-described by two springs and a damper [22,23], which can be described using the force balance equation as where σ and P 0 indicate the contact stress between pMMS and skin surface and barometric pressure, respectively. k a1 , k a2 , and c a are the elastic and damping components of the pMMS, respectively. From Eq. 2b and the derivative of Eq. 3 with respect to time, the dynamics between σ and P a can be described as Additionally, muscle-skin dynamics are known to have three distinguishable mechanical behaviors: the pure elastic, viscoelastic, and constant creep phase [24,25]. However, the constant creep phase only becomes dominant when deformation of the muscle-skin remains constant for long periods, which is difficult to be performed by voluntary muscle contraction. This high order property is further regressed by air pocket dynamics, thereby we assume that the constant creep phase is negligible. Then, the muscleskin dynamics are described by an elastic component, k m1 , serially connected to another elastic component, k m2 , and a damper, c m , to address the pure elastic and viscoelastic phases, or (k m1 +k m2 +c m s) and x m represents the η-dependent internal equilibrium point of the muscleskin dynamics. Eq. 5b was derived by the combination of Eq. 2b and 5a, putting Then, from Eq. 4b and 5b, the overall muscle-skin-pMMS dynamics are described as Because the air pocket is shielded by a double layer housing cell, the external perturbations to δP a , e.g. the muscular movements of surrounding muscles, and external forces, are negligible, i.e. δx 0 = 0. To enhance generality of the proposed model, the η-x m relationship is not restricted by the linearity assumption; the static nonlinearity g(.) is assumed between δx m and δη, or Then, as shown in Fig. 2a, the η-P a relationship is described as Therefore, δη can be estimated from δP a as Note that the inverse model of pMMS dynamics for the variable estimation can be described as the linear system G −1 sys and the static non-linearity property g −1 (.), or the Wiener system.

System identification
In this work, the system identification of pMMS is conducted by the revision of the Wiener system identification methodology addressed in [26]. Note that, because both G s and G m are first order dynamics with a single zero (Eq. 6a), G −1 sys (s) is described as a second order system with two zeros, or where a i and b j indicate coefficients. Then, the linear dynamics G −1 sys (s) and static non-linearity g −1 (.) of Eq. 9 are identified by the following process: 1. Estimate the best fit for the linear dynamicsĜ −1 sys (s) from δP a and δη. 2. Estimate output of the linear dynamicsŷ from the convolution ofĜ −1 sys (s) with δP a . 3. Fit a high-order polynomial g −1 p (.) betweenŷ and δη.
Note that the identified static non-linearityĝ −1 is determined as g −1 l (g −1 p (.)). In addition, because δP a has strong correlations with the muscle fascicle length l m and joint angle θ, we assume a second order polynomial for g p (.) during iterative system identification procedure.

Model verification of pMMS
The pMMS estimation of muscle fascicle length and joint angle is evaluated by comparing fascicle length and joint angle measurements via the ultrasound sensor and joint encoder as shown in Fig. 2b. The proposed model is then verified by measuring muscle motion of the tibialis anterior (TA). Five subjects (S1-S5) participated in this study as shown in Table 1. All experiments were carried out with informed consent at the Massachusetts Institute of Technology (MIT), under the approval of the Committee on the Use of Humans as Experimental Subjects (COUHES). The pMMS was affixed to the skin overlying the muscle using a compression sleeve or elastic wrap depending on the calf circumference and comfort of each subject. The pMMS placement was optimized to maximize the pressure sensor outputs during dorsiflexion. For system identification and evaluation of the pMMS estimations of the joint angle θ and muscle fascicle length l m , the joint encoder and ultrasound sensor were utilized as shown in Fig. 2b. For all trials, the subjects were asked to wear the joint encoder on their ankle to have the measurements of θ correspond to the pMMS outputs. Due to practical constraints surrounding ultrasound sensor placement alongside the pMMS, δl m is calculated from δθ values based on a lookup table capturing θ and l m which were identified prior to the trials without the pMMS. In addition, the δl m was normalized by the full range of δl m to evaluate over different individuals, referred to as m . To validate the proposed pMMS model, patients were asked to dorsiflex their foot, activating only the TA, from neutral position to full dorsiflexion. Full range of motion (ROM) that involves full plantar flexion would have introduced external perturbations from surrounding muscles into the measurement which will be detailed in the later portion of this report. In the first experiment, half-amplitude sweep trials, subjects were asked to dorsiflex at approximately 30, 60, and 90 percent of their full range of dorsiflexion at a steady pace. However, no specific joint trajectory was provided for the subjects to allow variance in the movements between each cycle. This range of movement enabled us to fit system parameters across a larger range of muscle movements. For the second experiment, half-frequency sweep trials, the subjects repeatedly performed dorsiflexion in the following pattern of paces: slow, fast, and hold; in which they stopped at full dorsiflexion for at least 3 seconds. Similar to the halfamplitude sweep trials, no specific trajectory guidelines were given to the subjects. Two data sets of each experimental scheme were collected. The system identification was conducted based on an integrated data set which consists one of two data set of both half-amplitude and half-frequency sweep trials. Model verification was conducted by comparing the estimates of the target variables, based on the identified system model, to their measured values of the trials that were not used in the system identification procedures.

External perturbations from surrounding muscles
In the previous experimental scheme, the pMMS model was evaluated under a limited ROM (half-amplitude and half-frequency sweep trials) to minimize external perturbations and focus on the sensor model verification. However, in practice, the muscular movements of surrounding muscles, e.g. peroneus longus (PL), tibialis posterior (TP), and antagonist muscles, e.g. gastrocnemius (GAS), may not be negligible and cause large estimation errors particularly at the extreme ends of ROM. A typical pMMS readings during full dorsiflexion and plantar flexion is demonstrated in Fig. 3a. During extreme plantar flexion, the pressure values do not uniformly correlate with the muscular movements due to external perturbations, which are marked as asterisks ( * ). For these values, because the pMMS cannot recognize whether the current pressure readings correspond to a small dorsiflexion or extreme plantar flexion, the external perturbations could cause significant errors in the both system identification and variable estimation as shown in Fig. 3c-f.
To investigate the effect of these external perturbations, all five subjects (S1-S5) were asked to perform two types of full range dorsiflexion and plantar flexion: full-amplitude sweep and full-frequency sweep trials. During the full-amplitude sweep trials, the subjects Fig. 3 External perturbations to the pMMS at extreme range of motion (ROM). The pMMS outputs do not uniformly correlate with the muscular movements at extreme ROM (circled with *). These perturbations are compensated by a differential-pMMS. a Typical single-pMMS outputs from the TA during dorsiflexion and plantar flexion, b schematic diagram of the differential-pMMS, c representatives of singe-pMMS (top) and differential-pMMS (bottom) outputs to varying joint angles, d representatives of singe-pMMS (top) and differential-pMMS (bottom) outputs to varying muscle fascicle length e representative linear plot of singe-pMMS (top) and differential-pMMS (bottom) estimation of varying joint angles, f representative linear plot of singe-pMMS (top) and differential-pMMS (bottom) estimation of varying muscle fascicle length repeatedly performed dorsiflexion and plantar flexion, gradually increasing the displacement from neutral position with each cycle. For the full-frequency sweep trials, the subjects were asked to repeat the same task as the half frequency sweep trial but with their full ROM (full dorsiflexion to full plantar flexion). No specific joint trajectory was given for any trials to allow their dorsiflexion and plantar flexion to be varied. At least two trials for each type of task were obtained for system identification and validation. The system identification and data processing were conducted in a similar manner to the previously described method.

External perturbation compensation by differential-pMMS
To compensate external perturbations from surrounding muscles, we use an additional pMMS as a reference value, and utilize the differential pressure of two pMMSs (differential-pMMS) to amplify the pressure changes by the muscular movement of interest and cancel out external perturbations as shown in Fig. 3b. This provides uniform positive or negative correlation with the target variables as shown in Fig. 3c-f. Thus, for the muscles that are affected by non-negligible external perturbations, the differential pressure δP diff can be used as to estimate δη where δP sen and δP ref indicate the pressure values of the pMMS on the muscle belly and reference point, respectively.

Identified system dynamics
The identified linear system dynamics, G −1 sys (s), are given in Table 2. See Supplemental information for subjectspecific dynamics. Large variances were found in the parameters of identified system dynamics between subjects. These variances are due to the sensitivity of the pMMS to sensor placements and tension of the sleeve to hold the sensor, which vary K p , k m1 , k m2 , c m as well as static nonlinearity between the target variable, δη, and the measured pressure, δP a .

Model verification of pMMS
The model verification results are shown in Fig. 4 and given in Table 3. For the joint angle estimation results, normalized root mean square (NRMS) errors were calculated by normalizing the root mean square (RMS) error in the ROM, which is referred to as e J . Meanwhile, because m is already the normalized variable, the RMS errors of m were used for the evaluation, which is referred to as e M . Additionally, the correlation coefficient (R-values) for the joint angle R J and muscle fascicle length estimations R M are calculated. The e J and e M are shown as 4.8 ± 0.5% and 5.7 ± 1.8% over the five subjects, respectively. Similarly, the R J and R M indicate strong linearity. The R J and R M are 0.990 ± 0.002 and 0.980 ± 0.014, respectively. Therefore, while the proposed model may simplify complex muscle-skin-pMMS dynamics, the current pMMS model was able to provide muscle fascicle length and joint angle estimation within 5-6% NRMS errors, suggesting validity of the proposed model for the motor tasks at different amplitudes and speeds tested here.

External perturbations from surrounding muscles
The pMMS performance under the external perturbations is shown in Fig. 5 and exhibited in Table 3. Comparing the results of half-amplitude and half-frequency sweep trials, the e J and e M of the full-amplitude and fullfrequency sweep trials drastically increased to 13.9 ± 1.7% and 11.9 ± 2.0%, respectively. Additionally, the R J and R M decreased to 0.884 ± 0.036 and 0.923 ± 0.016, respectively. As described in the "Methods" section, these deteriorations of pMMS performance are due to external perturbations from surrounding muscles as well as insensitivity of pMMS outputs to muscular movements around neutral position as shown in Fig. 3.

External perturbation compensation by differential-pMMS
The differential-pMMS system was verified by the fullamplitude and frequency sweep trials which demonstrates that the differential-pMMS effectively compensates for the external perturbations (Figs. 3 and 5, Table 3) and improves e J and e M to 7.0 ± 2.5% and 6.4 ± 1.2%, respectively. Also, the differential-pMMS estimation results show strong linearity to the target variables which R J and R M are 0.971 ± 0.018 and 0.979 ± 0.010, respectively. To further address the affect of differential-pMMS, one-way ANOVA with Tukey-Kramer multiple comparison tests were applied to e J , e M , R J , and R M of different pMMS configurations which are shown in Fig. 6. No significant differences were found between single-pMMS halfamplitude and half-frequency (S.H.H.) and differential-pMMS full-amplitude and full-frequency (D.F.F.) trials at all matrices including e J , e M , R J , and R M . However, significant differences were found between single-pMMS full-amplitude and full-frequency (S.F.F.) and the other two configurations at all matrices. These results indicate that the differential-pMMS configuration is effective in compensating for external perturbations and successfully provides reliable joint angle and muscle fascicle length estimation during the full ROM at the ankle.

Discussion
Throughout this work, we have investigated the characteristics and capability of the pMMS to provide muscle fascicle length and joint angle estimation. Using silicon tubing (Fig. 2b), the current pMMS design allows for implementation of the electronics outside of an fMRI scanner room, while keeping the air pocket attached to the patient, demonstrating its viability as a MRI-compatible sensor platform. Moreover, because the pMMS can be easily applied to the muscle of interest without being invasively attached to skin, the pMMS provides superior advantages in terms of practicality. Although highorder muscle-skin interaction dynamics were simplified by neglecting 'creep' dynamics, a single pMMS successfully provided comparable estimates to actual measurements of joint angle and muscle fascicle length during Fig. 4 Model validation results. a Representative joint angle estimation result from half-amplitude (top) and half-frequency (bottom) sweep trials (S1), b representative linearity plot of joint angle estimation (S1), c summary of NRMS estimation errors (error bars) and R-values (*) of joint angle estimation (S1-S5), d representative muscle fascicle length estimation result from half-amplitude (top) and half-frequency (bottom) sweep trials (S1), e representative linearity plot of muscle fascicle length estimation (S1), f summary of NRMS estimation errors (error bars) and R-values (*) of muscle fascicle length estimation (S1-S5). The dotted lines and solid lines in the estimation results indicate the measured values and estimation values from pMMS outputs, respectively. For the linearity plots, the solid lines and dotted lines indicate the mean and 1 S.D. of the pMMS linearity, respectively active dorsiflexion of all five participants within different speeds and amplitudes as shown in Fig. 4. To overcome large estimation errors resulting from perturbation of surrounding muscles during plantar flexion, shown in Fig. 5, we explored the integration of multiple pMMS to provide estimates of muscle fascicle length and joint angle. By having one pMMS at both the proximal and distal ends of the TA, we found that the estimated joint angle and muscle fascicle length from differential-pMMS readings were comparable to that of the measurements from the joint encoder and ultrasound imaging, expanding pMMS capacity over a full ROM. Therefore, the pMMS may lead to further understanding of CNS and muscle sensory feedback interactions, as well as enhanced rehabilitative strategies for populations with sensory-motor dysfunction.
Although the current pMMS poses sufficiency of its use in fMRI settings for real-time monitoring of a subject's task performance, more general application usage, such as human motor intent estimation for assistive device control and clinical measurements in rehabilitation Comparison of single-pMMS and differential-pMMS performances for full-amplitude and full-frequency trials. a Representative joint angle estimation result from full-amplitude (top) and full-frequency (bottom) sweep trials (S1), b representative linearity plot of joint angle estimation (S1), c summary of NRMS estimation errors (error bars) and R-values (*) of joint angle estimation (S1-S5), d representative muscle fascicle length estimation result from full-amplitude (top) and full-frequency (bottom) sweep trials (S1), e representative linearity plot of muscle fascicle length estimation (S1), f summary of NRMS estimation errors (error bars) and R-values (*) of muscle fascicle length estimation. The dotted lines and solid lines in the estimation results indicate the measured values and estimation values from pMMS outputs, respectively. For the linearity plots, the solid lines and dotted lines indicate the mean and 1 S.D. of the pMMS linearity, respectively settings, needs further investigation. Because muscle mechanical properties are also affected by the level of muscle activation, the current pMMS systems may be limited in providing estimates during isometric contraction or passive stretching movements which is commonly performed in rehabilitation and physical therapy settings. To overcome these potential limitations, the incorporation of EMG data in conjunction with pressure data may be used to enhance the usability of the pMMS outside of fMRI applications. In this scenario, EMG data can be utilized to decouple muscle fascicle length and muscle activity from pressure measures, applying advanced sensor fusion methodologies, e.g. Kalman filter [27], fuzzy logic [28].

Conclusion
In this work, the pMMS design and associated methodology are proposed. Through several diverse trials, investigating the pMMS's capability to estimate muscle fascicle dynamics under various movement frequencies and joint angle amplitudes, we report high prediction accuracy of joint position and muscle fascicle length. The pMMS poses as a reliable, MRI-safe, wearable sensor for measuring muscle fascicle length and dynamics in populations including persons with amputation or motor dysfunction because of its ability to perform independently from joint angle measurements. The application of the pMMS system in conjunction with fMRI data of brain activity Asterisks above bars indicates significant differences between different pMMS configurations. *p < 1 × 10 −3 , **p < 5 × 10 −4 , ***p < 1 × 10 −4 , ****p < 5 × 10 −5 , one-way ANOVA with Tukey-Kramer multiple comparison tests. Where no significance was seen, a p value is shown may lead to novel discoveries for the augmentation of rehabilitation and mobility. Future work will include the investigation of pMMS usage in non-MRI settings and the incorporation of EMG for more accurate muscle dynamic measures. Additionally, the pMMS will be tested during rehabilitation appropriate tasks such as passive stretching and weight-bearing activity.