Combining Ultrafast Ultrasound and High-Density EMG to Assess Local Electromechanical Muscle Dynamics: A Feasibility Study

Skeletal muscles generate force, enabling movement through a series of fast electro-mechanical activations coordinated by the central nervous system. Understanding the underlying mechanism of such fast muscle dynamics is essential in neuromuscular diagnostics, rehabilitation medicine and sports biomechanics. The unique combination of electromyography (EMG) and ultrafast ultrasound imaging (UUI) provides valuable insights into both electrical and mechanical activity of muscle fibers simultaneously, the excitation-contraction (E-C) coupling. In this feasibility study we propose a novel non-invasive method to simultaneously track the propagation of both electrical and mechanical waves in muscles using high-density electromyography and ultrafast ultrasound imaging (5000 fps). Mechanical waves were extracted from the data through an axial tissue velocity estimator based on one-lag autocorrelation. The E-C coupling in electrically evoked twitch contractions of the Biceps Brachii in healthy participants could successfully be tracked. The excitation wave (i.e. action potential) had a velocity of 3.9±0.5ms−1 and the subsequent mechanical (i.e. contraction) wave had a velocity of 3.5±0.9ms−1. The experiment showed evidence that contracting sarcomeres that were already activated by the action potential (AP) pull on sarcomeres that were not yet reached by the AP, which was corroborated by simulated contractions of a newly developed multisegmental muscle fiber model, consisting of 500 sarcomeres in series. In conclusion, our method can track the electromechanical muscle dynamics with high spatio-temporal resolution. Ultimately, characterizing E-C coupling in patients with neuromuscular diseases (e.g. Duchenne or Becker muscular dystrophy) may assess contraction efficiency, monitor the progression of the disease, and determine the efficacy of new treatment options.


I. INTRODUCTION
Muscles generate force through a series of fast electrochemical and electromechanical processes, enabling movement [1]. Muscle contraction is initiated by spinal motor neurons, which generate action potentials (APs). An AP travels through the motor nerve to the neuromuscular junction (NMJ, also referred to as the motor end-plate). Upon arrival The associate editor coordinating the review of this manuscript and approving it for publication was Yizhang Jiang . at the NMJ, the AP depolarizes the cell membranes of the innervated muscle fibers. The depolarization triggers a chain of cellular processes, which results in propagation of the AP along the muscle fiber in both directions towards the tendons with a velocity of 3-5 m s −1 [2].
A fiber consists of multiple myofibrils in parallel, and myofibrils consist of multiple sarcomeres in series. The sarcomeres are the smallest contractile units of a muscle, and following AP arrival, calcium will be released, allowing cross-bridges to form between actin and myosin filaments in VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ a sarcomere [1]. As the AP propagates longitudinally along the fiber, more sarcomeres will be activated and start to form cross-bridges. The cross-bridges allow generation of a contractile force [3]. This sequence of electrical and mechanical reactions is referred to as the excitation-contraction (E-C) coupling [4], and happens locally in the muscle fiber at a timescale in the order of milliseconds. The contractile force of a sarcomere is transmitted to the tendons through various pathways [5]. A clear distinction can be made between force transmission in the intracellular structures and in the extracellular matrix (ECM) [6]. Several studies indicate that muscle pathologies affect the force transmission and thereby the capacity to generate muscle force, i.e. the contraction efficiency [7]- [9].
Currently, there are no non-invasive techniques to assess muscle contraction and to track the E-C propagation. Several studies proposed to assess the electromechanical delay (EMD) as a measure for contraction efficiency. The EMD is defined as the time lag between onset of muscle electrical activity and force production at whole muscle scale (see Figure 1) [10]. However, the EMD at whole muscle scale provides no insight in the muscle fiber fast E-C dynamics, and thus the mechanisms of muscle contraction and force transmission.
Nordez et al. [11] proposed to quantify EMD with ultrafast ultrasound imaging (UUI) by assessment of muscle fiber and tendon motion onset in electrostimulated twitch contractions [11]. However, the velocities were averaged over the entire ultrasound field of view (FOV) of the muscle belly [11]- [15]. As a result, the E-C propagation was not tracked and the method provided only limited insight in the onset of muscle fiber contraction and fast dynamics of muscle fiber electromechanics.
In the current study, we propose a method that combines high-density surface electromyography (HD-sEMG) and UUI to follow the E-C propagation by simultaneously tracking the propagating excitation wave (the AP), and the subsequent local muscle fiber deformation caused by the contraction (i.e. a mechanical wave).
To aid interpretation of the experimental results, a novel nonlinear micromechanical muscle fiber model was required. The majority of muscle modelling studies either reduce the whole muscle to a single contractile element [16], [17] or represent the muscle as a multisegmental system comprising multiple Hill models in series and assume instantaneous muscle activation [18]- [21], thus neglecting E-C propagation. Other studies that modelled a whole muscle as a series of Hill models considered a time scale of seconds, limiting their use to study the E-C dynamics [22], [23]. Therefore, a multisegmental fiber model conceptually similar to the model by Morgan et al. was implemented [20]. By sequentially activating the individual sarcomeres, the fiber response during the first milliseconds of muscle contraction was simulated.
This paper describes and evaluates a novel experimental method to track the E-C propagation in vivo by . The classical EMD comprises the delay between stimulus and force production onset. Nordez et al. proposed a more complete characterization of the EMD by including motion onset of muscle fibers and the myotendinous junction (tendon-muscle fiber attachment point) [11]. Biceps motion onset and force onset timing adapted from [11]. There are currently no methods to assess the local electromechanical processes within the gray dashed square in vivo. The sequence of electromechanical processes is schematically represented in (b). After a stimulus, an Action Potential (AP) propagates along the fiber, causing the local fiber segments, i.e. sarcomeres, to be activated. Upon activation, cross-bridges will form within a sarcomere, causing the sarcomere to contract. This causes a local active motion onset in the fiber segment. The shown Excitation-Contraction (E-C) coupling timings are adapted from Sandow [4], and were based on single-fiber experiments. There is limited in vivo data available on the dynamics of the processes included in the brace under the arrow, indicated by the question mark.
combining UUI and HD-sEMG. In addition, a newly developed micromechanical muscle model is described and simulation results of the onset of contraction are presented. Finally, the similarities between the measured and the simulated E-C dynamics are discussed, and the implications of our observations are discussed.

II. METHODS
Three healthy men participated (22-25 years). The participants gave informed consent. The experiment was approved by the TU-Delft ethics committee (approval number 842). Data acquisition and processing was done in Matlab R2018b (The MathWorks, Natick, MA, USA).

A. INSTRUMENTATION 1) SURFACE ELECTROMYOGRAPHY
The action potentials were tracked using surface electromyography (EMG). An 8 × 8 square HD-sEMG electrode grid (inter-electrode distance 4 mm) was placed proximal with respect to the motor point over the short head of the Biceps Brachii (BB; see Figure 2b and d). A reference Ag-AgCl electrode was placed on the wrist. Unipolar signals were sampled at 2048 Hz using a 136-channel Refa amplifier (TMSi, Oldenzaal, the Netherlands).

2) ULTRASONOGRAPHY
Muscle contractions were imaged using a linear transducer array (L12-5 50 mm 256 element probe, pitch 0.2 mm; Philips, Bothell, WA, USA) connected to a Vantage-256 The solid black arrows represent data connections, the dashed lines represent TTL signal connections for synchronization of the equipment. The Verasonics ultrasound scanner sent a TTL trigger to the arbitrary waveform generator (AWG), the AWG sent a trigger (delay 20 ms) to the TMSI EMG amplifier and the muscle stimulator. (b) Zoom-in of the positioning of the measurement equipment on the Biceps Brachii (BB). The stimulator (blue rectangle) was positioned over the motor point, which is located on the distal part of the BB. The HD-EMG electrode grid (orange square) was placed proximal with respect to the stimulator and the ultrasound probe (red rectangle) was positioned over the center of the HD-EMG grid. (c) Ultrasound image of the Biceps Brachii (BB), the EMG electrodes are depicted on top of the ultrasound image. The red line denotes one of the M-modes that were used in the analysis. The acoustic shadow caused by the HD-EMG grid is visible directly beneath the HD-EMG grid. (d) Stimulator positioned on the motor point short head of the BB, the cathode was placed on the proximal side. The stimulator anode was positioned distal with respect to the cathode, in the same arm strap at a distance of 20 mm. (e) Ultrasound transducer placed on top of the HD-EMG electrodes. The arrows denote the coordinate system. research ultrasound scanner (Verasonics Inc., Kirkland, WA, USA). The scanner was set to transmit plane waves with a center frequency of 7.8 MHz, and sampled the raw radio frequency (RF) signals at 4 times the center frequency. Since the scanner consists of 128 channels per transducer connector, one image was formed using three overlapping synthetic apertures of 128 elements. Imaging depth was set to 16 mm, resulting in an effective frame rate of 5000 frames per second. Per contraction, 500 frames were recorded.

3) ELECTRICAL STIMULATION
Muscle contractions of the BB were evoked by bipolar electrical stimulation. Electrical stimulation results in the propagation of a compound motor action potential (CMAP), termed the M-wave. The stimulator (Micromed Energy, Micromed S.p.A., Mogliano Veneto, Italy) gave a single pulse with only negative phase, with a pulse duration of 200 µs and an amplitude between 5-10 mA. The stimulator cathode was placed over the motor point of the short head of the BB, which is located on the distal part of the BB [24], [25]. A graphical user interface (GUI) to monitor whether the stimulator did evoke a contraction was implemented in Matlab. Upon detection of a trigger, the EMG signal traces of both a distal and proximal EMG electrode pair were shown in the GUI window.

4) TRIGGERING AND SYNCHRONIZATION
To trigger the stimulator and to allow synchronization of the recordings of the ultrasound scanner and EMG amplifier in data processing, a transistor-transistor logic (TTL) pulse was sent from the ultrasound scanner to an Arbitrary Waveform Generator (AWG; DG1022A, Rigol, Beijing, China) at the start of each ultrasound acquisition. Upon receipt of the TTL pulse, the AWG transmitted a TTL pulse to the inputs of the EMG amplifier and stimulator after a delay of 20.00 ms. The delay allows the ultrasound scanner to start recording ultrasound frames prior to the contraction in order to have sufficient baseline without motion (see overview in Figure 2a).

B. PROTOCOL
Participants were seated in a comfortable, height-adjustable chair and were asked to rest their arm on a table in front of them (see Figure 2e). First, the motor point of the BB was detected by placing the stimulator on the distal portion of the BB, around (2/3) of the length of the BB. The stimulator was moved over the BB while repeatedly stimulating, until a point with notable activation was found (verified by the amplitude of the EMG signal traces in the GUI). When found, the stimulator was fixated using an elastic arm strap. The stimulator was set to a low current, such that only muscle fibers directly underneath the cathode were contracting, and the contraction force did not overcome joint friction. Before placement of the EMG electrodes, the skin on the upper arm and wrist was shaved and cleaned with an abrasive skin preparation gel and alcohol to reduce skin impedance. Next, the ultrasound transducer was placed on top of the HD-sEMG VOLUME 9, 2021 electrode grid, and fixated using a holder (MG70003, Noga Engineering & Technology Ltd., Shlomi, Israel). The transducer did not compress the muscle, and made contact with the HD-sEMG grid through acoustic transmission gel. The transducer position was adjusted such that the full length of the HD-sEMG grid was in the field of view (FOV) and the majority of the muscle fibers in the FOV were aligned with the longitudinal axis (verified by conventional real-time B-Mode imaging).
The experiment consisted of three sets of ten electrical stimulations, with an inter-stimulus interval of around 5 s and a 5 minute break between sets. After each break, the stimulator position was checked to still elicit muscle activation and adjusted if needed.

C. PROCESSING 1) MUSCLE FIBER CONDUCTION VELOCITY ESTIMATION
EMG signals were first spatially filtered using the longitudinal single differential (longitudinal direction, Figure 2c and e) to obtain bipolar signals [26]. Next, the signal was band-pass filtered using a bidirectional 3rd-order Butterworth filter (cutoff frequencies 10-250 Hz).
The measured bipolar signal per electrode over time is given by y ck (t), with k = 1 . . . K the electrode row (proximal to distal) and c = 1 . . . N c the electrode column (medial to lateral).
To find the average muscle fiber conduction velocity (MFCV) of the CMAP, a multichannel maximum likelihood estimator in the frequency domain was used [27]. This resulted in an estimate of the propagation delayθ per recording. Then the velocity of the CMAP, V cmap , was found by V cmap = d/θ, with d the inter-electrode distance.
The algorithm does not take into account misalignment of the electrode grid with respect to the muscle fibers. However, the Biceps Brachii is a fusiform muscle, and as seen by B-Mode imaging, the misalignment was small. For small angles of misalignment (i.e. < 5 • ), the influence of misalignment on the estimated MFCV was considered negligible.
The negative peak in the CMAP was considered as the depolarizing peak in the EMG. The arrival time of the CMAP T arr-cmap was found by minimizing the following equation, To not be limited by the sampling frequency, linear interpolation of y ck (t) was performed during the minimization. The optimized variable t 0 resulted in the arrival time T arr-cmap (depicted in Figure 3).

2) ULTRASONOGRAPHY
Raw RF channel signals were processed to beamformed In-phase Quadrature data using conventional delay and sum beamforming, resulting in S(x, z, t i ). The EMG and ultrasound data were spatially aligned using the acoustic shadow present due to the electrode grid (see Figure 2c). Tissue velocities v z (x, z, t i ) along the ultrasound beam axis (i.e. z-axis, depth Figure 2c) were determined using one-lag autocorrelation speckle tracking [28]. Since only axial velocities are determined, the obtained velocities do not directly reflect muscle fiber longitudinal deformation, but rather muscle tissue thickening and thinning. The tissue velocities v z were temporally filtered using a bidirectional 3rd-order Butterworth low-pass filter (cutoff frequency 100 Hz). Next v z was differentiated to obtain axial tissue acceleration a z , with F s the frame rate. The mechanical wave was tracked using a normalized Radon transform of the acceleration a z [29]. To find the onset of tissue motion, only a z (x, z, t i ) within a time window of moment of first detection of EMG activity and the second (positive) peak in the EMG was used, Furthermore, only data in a region of interest (ROI), chosen directly below the EMG grid, was used to track the mechanical wave (horizontal lines in Figure 4). Three horizontal M-modes in the ultrasound image at depths where a bundle of muscle fibers was clearly visible were chosen manually for analysis. Acceleration data along these M-modes (averaged over ±0.5 mm in depth), within the ROI and the time of window of interest was used to track the contractile mechanical wave (ConW), resulting in a velocity V conw and time of arrival T arr-conw . Peaks in the Radon transform corresponding to a velocity outside V cmap ± 2 m s −1 were regarded out of physiological range and excluded.

D. FIBER MODEL
A multisegmental fiber model conceptually similar to the model by Morgan et al. was implemented [20]. The model is able to simulate the onset of muscle fiber contraction, and served to interpret the experimental results. A sensitivity analysis was carried out to see the effect of varying AP velocity on the mechanical fiber response. The full model description and the results of the sensitivity analysis can be found in the Supplementary Materials. The muscle fiber was modelled as a large chain of n = 500 sarcomeres in series. In the following, the subscript i denotes sarcomere number. The force in sarcomere i is described by the Hill muscle model, relating the sarcomere length i (t), velocity˙ i (t) and activation a i (t) to force F s,i (t). Only the parallel and contractile elements were considered for an individual sarcomere, separating the force in an active part due to cross-bridge interactions, and a passive part caused by connective tissue in the fiber. Hence the force per sarcomere is given by, where F PE denotes the passive parallel elastic sarcomere force, F CE the contractile element sarcomere force and f a the sarcomere activation function. The sarcomeres are separated by boundaries, i.e. Z-disks, which are represented by a 1D array of m = n+1 point masses (i.e. nodes), where boundary j has position x j (t). The first boundary is fixed, and the last boundary x m (t) is connected to a tendon, modelled as a linear spring. The tendon corresponds to the series element in the Hill model. The position of the other end of the tendon is given by x t , and is fixed during muscle contraction.
To model the propagation of the action potential, the sarcomere activation a i (t) is dependent on the initial sarcomere position and the velocity of the action potential V ap . The sequential activation of the sarcomeres is modelled by a linear increasing activation function f a,i . The Hill model uses an input signal in the interval [0, 1], where f a,i = 0 corresponds to an inactive sarcomere and f a,i = 1 to a fully activated sarcomere. Fiber relaxation was not simulated. The activation function of sarcomere i is given by, with t d,i the delay for activation onset and t a,i the time required to reach full activation. The delay for activation onset t d,i for sarcomere i due to propagation time for the action potential along the fiber is given by, with optim,i the optimum length of the sarcomere, t ss a general delay term and V ap the velocity of the action potential. A damping force was implemented on the nodes between the sarcomeres, to avoid expected high frequency oscillations due to the low fiber mass. The damping force on node x j (t) is given by, with c damp,j the damping coefficient for node j. Since the high-numbered boundaries (close to the tendon) are expected to have higher velocities than the low-numbered boundaries (the first boundary is fixed), the damping coefficient for boundary j is scaled by its initial position x 0,j , Experimental results for three participants, M denotes mean, SD standard deviation, SR success rate. The velocity of the action potential was found for all trials of participant 1 and 3, and in 20/30 for participant 2 (CV found). T arr-cmap denotes the arrival time of the depolarizing peak in the CMAP. The number of successful fits and the results for Radon transforming acceleration data to find the velocity V conw of the contractile wave are presented for three M-modes at various depths. T arr-conw denotes the arrival time of the contractile wave. Note that T arr-conw is sometimes similar to T arr-cmap and sometimes after T arr-cmap . This is caused by the Radon tracking method, which in some cases found the precontractile motion (blue area in Figure 4 and 5), and in other cases the later occurring higher amplitude motion. The difference in motion used for Radon fitting did not have a large influence on the detected contractile wave velocity V conw .
This scaling of the damping can be justified by considering that during a whole muscle contraction the intramuscular fluid and the surrounding elastic structure will move together with the fiber of interest. Hence, in reality the fluid will have a lower relative velocity with respect to the Z-disks then in the modelled fiber, where the surroundings are fixed.

SIMULATION PROCEDURE
Numerical integration of the fiber model was done using Matlab 2018b (The MathWorks, Natick, MA, USA). Since the system has very low mass, the stiff ode15s solver was used for numerical integration. Tolerances were set to 1 × 10 −10 and integration time was set to 20 ms. Since ODE solvers cannot handle discontinuous functions well, continuous functions were fit to the piecewise continuous functions from the Hill model (described in the Supplementary Materials). In the experiment, the elbow was near full extension, causing strain on the muscle fibers, and thus some passive tension in the muscle. To simulate the model in a similar scenario, the initial condition for the fiber was set to zero velocity and an initial strain ε 0 of 7.5 %, causing passive pre-tension in the fiber.

A. EXPERIMENTAL RESULTS
Transcutaneous electrostimulation of the motor point of the Biceps Brachii successfully evoked muscle contractions in all participants, verified by the EMG CMAP recordings. All measured CMAPs propagated from distal to proximal, away from the motor point. Conduction velocities were estimated successfully for 80/90 trials, see Table 1. The EMG measurement was highly repeatable, illustrated in Figure 3. The time of arrival of the depolarizing peak in the CMAP, T arr-cmap , was in all trials well after the time instant of stimulus.
Propagating mechanical waves were observed in the tissue velocity profiles, as well as in the tissue acceleration profiles (see Figure 4, 5 and S4). The Radon transform of the acceleration data within the ROI successfully tracked the mechanical waves in 83/90 trials. All waves propagated from distal to proximal. The tracking results varied with M-mode depth, and the total of 90 Radon transforms per participant resulted in a success rate of 57% for participant 1, 73% for participant 2 and 59% for participant 3 (see Table 1). Both the velocity, V conw and time of arrival, T arr-conw , of the mechanical wave had larger spread than V cmap and T arr-cmap , as can be seen by the standard deviations. In some cases, the arrival time of the mechanical wave preceded the arrival of the CMAP (e.g. participant 1, M-modes 10 and 14 mm, see Fig. S3 and S4). In most trials, the velocity and acceleration profiles showed a small positive peak first (indicating tissue moving deeper), which was followed by a negative peak (indication motion in superficial direction) depicted in Figure 5 and S8. The peaks in the acceleration profiles where more distinct than the peaks in the velocity profiles.

B. FIBER MODEL RESULTS
The fiber model was able to simulate the propagation of the action potential, and the subsequent contraction onset of individual sarcomeres. When looking at what occurs locally in the muscle fiber, we see that there is non-uniformity in sarcomere length and velocity (see Figure 6). The standard deviation in sarcomere length was larger for the high-numbered sarcomeres. From Figure 6 it becomes clear that the standard deviation in sarcomere velocities is high at the moment where the respective groups of sarcomeres start to be activated. Looking at the sarcomere lengths in Figure 6, it can be seen that the sarcomeres are first elongated, and later start to shorten when the sarcomeres are locally activated. This effect is strongest for the high-numbered sarcomeres.

A. EXPERIMENTAL RESULTS
The current study demonstrated that by simultaneously using high-density surface electromyography (HD-sEMG) and ultrafast ultrasound imaging (UUI), fast electromechanical muscle dynamics can be tracked. Propagating CMAPs were detected in all trials and propagated from distal to proximal. The conduction velocity V cmap could be estimated in 80/90 trials. During the last 10 trials of participant 10, The CMAP is propagating from left to right, i.e. distal to proximal. The orange triangle denotes the depolarizing peak in the CMAP. From 10 ms a positive acceleration wave propagates from left to right, likely caused by precontractile motion. The blue area is immediately followed by a red area, denoting higher amplitude negative acceleration, likely caused by the onset of muscle contraction. Around 24 ms the acceleration flips sign again, and a positive acceleration wave propagates from distal to proximal. Notice the asymmetric colorbar, to make the positive acceleration stand out more, showing the precontractile response more clearly.
the stimulator position was most likely slightly off, causing a low amplitude CMAP compared to the stimulation artefact. This biased the estimated V cmap to high conduction velocities outside of physiological range [30]. Still, mechanical waves could be detected in these trials. CMAP conduction velocities were found to be within the physiological range, as reported in literature [1], [31]. The highly repeatable EMG recordings were comparable to other studies that used transcutaneous muscle stimulation [31]- [34]. The stimulation was done on the motor point with low current (< 10mA), and consequently, only superficial muscle fibers were stimulated. This was visually verified, as only very small movements of the skin were observed. Similar results of superficial stimulation were reported by [28] (see ref. Fig. 5 and 13).
The tracked mechanical waves had a velocity in the same order of magnitude as the CMAP conduction velocity. The mechanical wave also propagated from distal to proximal, with a mean (SD) velocity of 3.52 (0.89) m s −1 . The mean and standard deviation of the mechanical wave velocities were similar to values reported in previous studies [28], [35], [36]. However, the referred studies did not look at axial acceleration and instead used the axial tissue velocity data to track the mechanical wave. Furthermore, the referred studies did not report the rapid sign change in the velocity and acceleration profiles in the first few milliseconds after stimulation like we observed. The low amplitude peaks in opposite direction in the velocity profile prior to the larger motion (i.e. the blue positive areas in the velocity panels of Figure 4 and Fig. S6-8), were also reported by [11] (see ref. Fig. 3). In the study of Nordez et al. [11], this initial peak was absent in the distal portion of the muscle (i.e. close to the innervation zone) but prominent in the proximal part, comparable to the result shown in Figure 4. The onset of large tissue motion was in all trials well after the onset of EMG activity, depicted in Figure 4. In [35] the large contraction motion was also succeeding the EMG activity, however in the referred study EMG was recorded in separate trials prior to the ultrasound measurements.

MECHANICAL WAVE INTERPRETATION
In the current analysis, the axial tissue velocity and acceleration were used to track mechanical waves. The axial velocities do not directly reflect longitudinal fiber deformation, but instead were interpreted as fiber thinning and thickening. Under the assumption that muscle fibers are incompressible [37], the shortening of a muscle fiber will increase its diameter. This complicates interpretation of the structural changes happening locally in the fiber during the contraction, as we are looking at a secondary effect. However, Lopata et al. found a high level of similarity in strain across longitudinal and axial directions during muscle contractions recorded at low frame rate with a 3D ultrasound transducer, indicating VOLUME 9, 2021 that axial deformations can be used as a measure for fiber elongation [38]. It must be noted that the high similarity in strain in different directions may not be true for pennate muscles. Moreover, the mechanical waves were detected with a velocity in the same range as the CMAP conduction velocities. For these reasons, we interpret the mechanical waves as contractile waves (ConW). The ConW is thought to depict the onset of local fiber motion, initiated by the arrival of the AP. This claim is further supported by the fact that the high amplitude motion occurs after AP arrival.
The CMAP and ConW arrival times (Table 1) do not reflect the full dynamics of the excitation-contraction (E-C) coupling, for two reasons. (1) The CMAP arrival was determined based on the spatio-temporal location of the depolarizing peak in the surface EMG recordings (denoted by the arrow in Figure 3). However, (surface) EMG measures potential differences caused by the AP [1], which will initiate the calcium dynamics, resulting in the formation of cross-bridges and thus force generation (i.e. the E-C coupling). Therefore, the actual moment of AP arrival occurs earlier than the depolarizing peak and the onset of local cross-bridge formation, i.e. the moment of AP arrival coincides with the onset of EMG activity. (2) The peak in axial acceleration does not correspond to the onset of local contraction. The moment of motion onset is actually given by the onset in acceleration, and not the peak. However, detecting the onset of acceleration was not robust due to spatial variations in amplitude (depicted in Fig. S6 and S7). This also explains the differences in T arr-conw in Table 1. For these reasons, temporal peak detection of CMAP and ConW do not completely characterize the E-C timing.

IN VIVO: LOCAL SARCOMERE PULLING
We hypothesize that this precontractile response is caused by earlier activated fiber segments closer to the innervation zone, which will start to contract and pull on neighboring fiber segments. This pulling causes the neighboring (still inactive) fiber segments to elongate and thin. The pulling effect was in some trials observed prior to CMAP arrival, further indicating that this local motion is likely not caused by the active contraction of fiber segments at that location. A similar precontractile response was observed by [11] (using UUI) and by [39] (using a single stethoscope microphone). In both studies, this phenomenon was ascribed to latency relaxation (LR). LR was first observed by Rauh in 1922, and is described as a brief period in which the fiber tension momentarily drops, causing a muscle fiber to elongate prior to contracting [4]. Since the AP propagates, we expect that the LR effect will propagate as well, resulting in the propagating precontractile response that we observed. LR has been measured extensively in isolated frog muscles, but recently evidence of a similar precontractile response in vivo was found using near-infrared techniques [34]. Using this optical technique a precontractile response temporally coinciding with the EMG activity was found. Since it is currently unknown how the ultrasound and optical response are linked, it remains elusive whether the precontractile response detected in the tissue acceleration is indeed caused by LR, and corresponds to the optical response measured by [34]. However, if LR causes a rapid drop in local fiber tension, this would result in a more clearly perceivable sarcomere pulling effect, supporting our hypothesis.

EXPERIMENT METHODOLOGICAL CONSIDERATIONS
Although longitudinal tissue velocities can be estimated, as shown by [35], [38], the spatial resolution in longitudinal direction is lower than the resolution in axial direction [28]. RF-based speckle tracking in axial direction results in a sub-wavelength spatial resolution, whereas in longitudinal direction the resolution is dependent on the ultrasound element distance (typically in the order of one wavelength). Furthermore, there is no phase information available in the longitudinal direction. To detect the onset of fiber motion during the initial phase of contraction, the high temporal resolution provided by UUI is essential. There have been various studies combining conventional ultrasound imaging and HD-sEMG to study muscle contractions (e.g. [40]- [43]), however, our study has an unprecedented temporal resolution in simultaneous recording of EMG and ultrasound, caused by the use of UUI. It is expected that our method enables new research into the E-C coupling in vivo.
For participant 1, the arrival time of the mechanical wave preceded the arrival time of the CMAP. This does not conflict with the presumption that these mechanical waves represent contractile waves, since this difference in arrival time seemingly resulted from the Radon tracking method. The wave tracking method could fit on either the first occurring peak (the presumed precontractile response), or on the higher amplitude, later occurring motion. In some trials this second motion was not a distinct peak, and therefore the Radon tracking method only detected the first peak. Further improvements in ultrasound data processing should be made to interpret the mechanical waves and to understand the physiological origin of the waves present during the onset of muscle contraction.
Comparing the velocities of the CMAP and ConW in Table 1, the waves' velocities are not the same. The mechanical wave velocity also has a larger standard deviation. There are two main reasons for these differences. (1) Presumably, the mechanical waves are composed of overlapping contractile waves and shear waves. The shear waves are induced by the contraction, and propagate through the stiff connective tissue (i.e. ECM) in the muscle. Since this tissue is stiffer than the contractile elements, the shear waves propagate faster, resulting in a spatially distorted wavefield. (2) There is more uncertainty in the ultrasound measurement and processing than in the EMG measurement and processing. To obtain tissue motion estimates from ultrasound data, and find the velocity of the mechanical wave is more challenging than finding the velocity of the CMAP. One factor is that the image quality varies within the image, because the orientation of the fibers with respect to the probe differs along the probe aperture.

B. FIBER MODEL SIMULATION RESULTS
A micromechanical multisegmental fiber model was developed. The model allowed simulation of the local contractile response to a propagating AP, i.e. the propagation of the excitation-contraction (E-C) coupling. Considering the general trend of the fiber contraction, the results were visually similar to those reported by Morgan et al. [20] (see S2). Morgan et al. modelled the fiber as 400 sarcomeres in series, with different initial lengths, and activated all sarcomeres simultaneously. The time to reach maximum force differs substantially between our results and those reported by Morgan et al., which can be explained by the parameters used in the model. In our model, the parameters that influenced the fiber contraction behaviour most were varied in the sensitivity analysis, i.e. nodal damping and AP velocity (see Table S1).

SARCOMERE LENGTH DISPERSION
Comparing the simulation results to other studies using multisegmental models (albeit other studies used simultaneous sarcomere activation), we observed similar results in sarcomere length non-uniformities, i.e. dispersion [18]- [20]. In our study, the sarcomere length dispersions were observed for all AP velocities used in the sensitivity analysis. We observed a difference in the sarcomere length distribution between lowand high-numbered sarcomeres. Although this has not been shown explicitly in simulation studies, Moo et al. measured similar non-uniformities in sarcomere length (i.e. SL dispersions) in vivo during isometric mice tibialis-anterior muscle contractions [44]. It was found that SL dispersion can be up to 1.0 µm for 30 serial sarcomeres. In addition, it was found that during activation, the sarcomeres at the middle of the muscle shortened by 7%, while the sarcomeres at the distal end shortened only 0-3%. Although these results are not directly comparable to the difference in low-and high-numbered sarcomere lengths in our simulations, it does highlight that sarcomere length varies locally throughout a fiber.

IN SILICO: LOCAL SARCOMERE PULLING
Upon closer examination of the sarcomere lengths during the initial phase of the contraction, we can see that some sarcomeres are first elongated prior to contracting (see Figure 6, S3 and S4). This happens in high-numbered sarcomeres and is caused by a pulling force of low-numbered sarcomeres that are activated earlier, similar to the local pulling effect observed in the experiment. A similar initial elongation of sarcomeres was seen by Stoecker et al. [19] (see ref. Fig. 3), though the authors did not explicitly address this phenomenon since their interest was on a longer time scale (up to multiple seconds).
A direct consequence of the local sarcomere elongation prior to contraction is the build-up of passive tension in the fiber. The force generated by actively contracting sarcomeres is first stored as elastic energy in the passive elastic connective tissue, and is later released in the form of kinetic energy. Another observation is that there is negative work done by sarcomeres that are elongated prior to contracting (see Table 2). As the elongating sarcomeres are first gaining eccentric velocity, this velocity first has to be reversed to start contracting. Due to the low inertia of the sarcomeres, this happens quickly. In the sensitivity analysis it became clear that the negative work was predominantly done by the passive tissue, however, for some parameter sets the active sarcomere force also contributed to the negative work done. This local negative power output was also present in the multisegmental model of van Leeuwen and Kier [22]. This effect can only be seen in multisegmental models, since these local sarcomere interactions cannot be studied in single contractile element models.
Another effect of this local pulling is that the sarcomeres are first being pulled more towards the descending limb of the force-length relation in the Hill model. This causes a decrease in the active force generating capacity of the sarcomere. When considering only the static case, this would indeed decrease the contraction efficiency of the fiber. However, in the dynamic case, the local eccentric velocity of the sarcomeres increases the force generating capacity following the force-velocity relation defined in the Hill model. As such, VOLUME 9, 2021 the Hill model makes it difficult to predict what exactly is happening locally in the fiber concerning contraction efficiency.

MODEL METHODOLOGICAL CONSIDERATIONS
The proposed model structure represents a single muscle fiber. Consequently, only series pathway force transmission can be studied, since fiber interactions could not be included. In the single fiber simulations, high contractile velocities and accelerations were observed, which are unlikely to occur in whole muscles. Increasing the density of the sarcomeres did limit the rapid velocity changes, but did not alter the general behaviour (i.e. the local pulling effect was still present). A whole muscle consists of fascicles, and fascicles consist of multiple fibers in parallel. These fibers have interconnections and are surrounded by connective tissue, resulting in various force transmission pathways [5]. A next step with this model would be to add multiple fibers in parallel, to further study the transmission of force during the onset of muscle contraction. This would also make the model suitable to study specific neuromuscular diseases, e.g. further our understanding of the absence of functional dystrophin in Duchenne muscular dystrophy.

C. IMPLICATIONS
The experiment and simulations both provided data that indicate local pulling of sarcomeres during the onset of contraction. It is currently uncertain whether the simulated data and the experimental data show the same phenomenon. The model simulated only a single fiber, hence no interactions with other fibers and connective tissue were considered. However, due to the used experimental conditions, it is not unlikely that the local pulling effect was present in vivo and in silico.
Latency relaxation (LR) in single muscle fibers has been shown extensively in previous studies [4]. The arm position used in the experiment causes some pre-stretch on the muscle fibers. Increased initial sarcomere length has been shown to amplify the tension drop in LR [45].
The similar behaviour in sarcomere pulling seen in silico cannot be originated by LR, since LR is not described in the Hill model. A future improvement of the model would be to make a more distinct separation in contributions of active muscle tissue and passive connective tissue. As such, a physiological model instead of a phenomenological model would be better suited to make such separation, at the cost of computational complexity [3], [18], [46].
The current study highlights that there are many unknowns in what happens locally in muscle fibers when the AP arrives, and how the fiber will respond locally to this excitation, i.e. the local E-C coupling. The presented novel method to image the E-C propagation seems promising to further investigate local fiber behaviour. The method needs to be standardized, and used to characterize the influence of various experimental conditions, such as the effect of joint position (i.e. effect of initial passive pre-tension and sarcomere length), the effect of stimulation current and the effect of measurement location with respect to the innervation zone. The in this paper presented results and the results of previous studies highlight that assuming a uniform sarcomere length and activation over the whole muscle is likely to result in erroneous prediction of muscle behaviour during contractions, especially during rapid transients [44]. Ultimately, the proposed experimental method can bridge the gap in fiber behaviour at the micro scale, and whole muscle behaviour at the macro scale [18].

V. CONCLUSION
• The novel experimental method to track electromechanical dynamics in skeletal muscle with high spatio-temporal resolution successfully detected propagating mechanical waves in response to propagating excitation waves (i.e. compound motor action potentials, M-waves).
• The newly developed multisegmental fiber model allowed to simulate the contractile response to a propagating excitation wave. Simulated fiber contractions showed qualitatively similar results as the experiment, but validation of multisegmental muscle models during the onset of contraction is warranted.
• Both the experiment and simulation resulted in precontractile motion that may be caused by earlier activated sarcomeres that pull on sarcomeres that were not yet excited. This effect was presumably amplified in the experiment by latency relaxation.
• Although many questions remain about the onset of muscle contraction, we showed a promising new technique that can provide additional insights for the integration of micro-scale and macro-scale measurements, and further our understanding of muscle behaviour. As a biologist with a Ph.D. degree in mechanical engineering, he (co-)developed experimental techniques to understand muscular and musculoskeletal function. In cooperation with the Delft University of Technology and MOOG this resulted in robotic wrist and ankle manipulators that were applied for experimental diagnostics of spasticity and from which a new approach of precision orthotics came forward. He is a Co-project leader of a stiffness compensating ankle foot orthosis is now being tested within the NWO-TTW Stiffness as a Needed-project.