Manual and semi-automatic determination of elbow angle-independent parameters for a model of the biceps brachii distal tendon based on ultrasonic imaging

Abstract

Tendons consist of passive soft tissue with non linear material properties. They play a key role in force transmission from muscle to skeletal structure. The properties of tendons have been extensively examined in vitro. In this work, a non linear model of the distal biceps brachii tendon was parameterized based on measurements of myotendinous junction displacements in vivo at different load forces and elbow angles. The myotendinous junction displacement was extracted from ultrasound B-mode images within an experimental setup which also allowed for the retrieval of the exerted load forces as well as the elbow joint angles. To quantify the myotendinous junction movement based on visual features from ultrasound images, a manual and an automatic method were developed. The performance of both methods was compared. By means of exemplary data from three subjects, reliable fits of the tendon model were achieved. Further, different aspects of the non linear tendon model generated in this way could be reconciled with individual experiments from literature.

Introduction

The prediction of human elbow joint movements can be achieved by means of biomechanical models of the skeleton, soft tissue, and muscles. Muscle forces are actively generated by muscle contraction, and transmitted to bones via tendons. The force generated by a muscle depends on its contraction velocity, length, and neural activation. The neural activation is usually measured indirectly in the form of surface electromyograms (sEMGs). The contraction velocity and length of the muscle can then be estimated based on biomechanical models which make use of sEMG signals. In case of a musculoskeletal model including tendons, the muscle length depends on the skeletal geometry and tendon length. The length of a tendon depends, among other, non linearly on the exerted force (elasticity) [1]. In several in vitro experiments the elastic properties of mammalian tendon tissue have been measured [1–3], but data on the distal biceps brachii tendon is sparse [4, 5]. Further, for in vivo elbow movement prediction, an in vivo method to access biomechanical tendon parameters is preferable, as in vitro experiments with tendons of living subjects requires interventions in the body. In vivo approaches to analyze the elastic tendon properties are either invasive by implanting sensors into the tissue or non-invasive by monitoring the tendon shortening by medical imaging and tendon force estimation based on external force measurements [6]. In this context, ultrasonography can be used as an imaging technique [3, 7–10]. For the upper arm flexor muscles there are currently only a few studies on in vivo tendon properties as opposed to those studies which focus on tendons of muscles in lower extremities used for walking [11]. The former upper arm studies focus on the distal biceps brachii tendon [9, 10] but with an emphasis on the analysis of linear stiffness parameters for medical applications.

As stated above a model, which predicts the elbow movement, depends on the tendon length. Therefore, the behavior of the biceps brachii tendon is of interest, as the biceps brachii is a key contributing muscle in elbow flexion [12, 13]. Therefore, in this work, the goal was to obtain a configurable model for the distal biceps brachii tendon which reflects the non linear stiffness material properties of the tendon and can be used in conjunction with a muscle model for model-based elbow movement prediction.

In order to obtain stiffness model parameters for living subjects ideally the tendon length over tendon force would be measured. Here, the displacement of the distal biceps brachii myotendinous junction parallel to the humerus was observed by means of ultrasonic imaging at different loads, for different effective lever arms (due to varying joint angles) and while the subject maintained a fully supinated hand. As the tendon load could not be measured directly, it had to be estimated based on a skeletal model. An abstract skeletal model including a simplified biceps brachii tendon path was used to predict the tendon force corresponding to the force measured at the wrist for each elbow angle. The tendon displacement parallel to the humerus was extracted from the ultrasonography images. For this purpose, a manual and an algorithmic method were used. The algorithmic method was compared with the manual method for data extraction. It was found that the presented algorithmic method was partially not applicable resulting in an semi-automatic data extraction. However, the semi-automatic extraction greatly reduced the effort of feature extraction.

In the following, a general overview on mammalian tendon properties is given before a non linear tendon model is introduced. The presented stiffness model was fitted onto the manually and (semi-) automatically extracted tendon displacement dependent on the predicted tendon force.

Mammalian tendon modeling

In general, tendons connect skeletal muscles to bones. Tendons insert into a bone via the enthesis and are connected to the muscle via the aponeurosis. The whole junction complex between muscle and aponeurosis is called myotendionous junction. In case of the distal biceps brachii myotendinous junction the aponeurosis additionally connects the biceps brachii with the forearm fascia [14].

The role of a tendon in joint dynamics is determined by the geometry and material parameters [15]. Independent of the role, each tendon is understood to have elastic, movement damping properties [16]. Additionally creepage [17] and friction with the surroundings [18] were observed. A tendon model covering the aforementioned properties (except creepage) is shown in Fig 1B. The model in Fig 1B is composed of principal mechanical units.

[Figure omitted. See PDF.]

(A) Pictogram of the biceps brachii (in red) and its lower tendon (in orange) attached to the radius. (B) A general tendon model including the tendon weight m_T, non linear stiffness F_k(Δl), non linear damping , and friction to surrounding tissue . (C) This model is obtained by applying static conditions to the model shown in (B).

If static conditions where assumed, motion-dependent, i.e. also time-dependent, properties would have no influence on the behavior of the presented model. If the assumption of static conditions is applied the tendon model depicted in Fig 1B is reduced to the elastic tendon properties and friction, as (Fig 1C). Furthermore in this work, friction is assumed to be negligible F_fric ≈ 0. As a consequence, only the elastic property of the tendon is further discussed. Elasticity is described as a force resisting changes in length. This is modeled by a function F(Δl). Where(1)is the difference between the actual length l and the slack length l₀. Usually, elastic properties are expressed as force over strain where strain is defined as . This implies knowledge on the slack length l₀. The slack length of a tendon in in vivo conditions is not necessarily known. Depending on the tendon in question and the experimental equipment that is available, it is likely that the slack length l₀ has to be guessed. This work presents a method for myotendinous junction displacement estimation in vivo in which Δl is measured. Therefore, the tendon has not to be modeled based on its resting (slack) length which is normally not available in vivo.

Nonlinear static mammalian tendon properties

Mammalian tendons, which consist of collagen structures, barely resist shortening therefore,(2)is assumed. The tendon resists lengthening in a non linear behavior, as demonstrated in Fig 2A. This behavior is due to the collagen structure in which the collagen strains do engage at different lengths [15]. When the collagen structure is fully engaged the elastic behavior becomes linear. The tendon length at this state is called Δl_toe. The region 0 ≤ Δl_T < = Δl_toe is referred to as toe region as defined by [1]. If the tendon is at Δl_T = 0, it is at its slack length, the length at which there is virtually no resistance to lengthening.(3)

[Figure omitted. See PDF.]

In (A) the tendon force is shown over the change in tendon length in orange. As described in Eq (10). The linear force length relationship is extended with a dashed line. In (B) the stiffness k_T of the exponential tendon model is shown over the tendon force . The tendon model parameters are marked in the plot. The non linear stiffness is configured by (c₁ and c₂). As shown in Eq (6) c₁ is the change in stiffness k_T over a unit change in F. the parameter c₂ is the stiffness at . At the force F_toe, the model continuously transitions into a linear stiffness function.

For modeling the tendon force in relation to changes in length there are two major models. The first is a polynomial tendon model as described by Buchanan et al. [19]. The second is an exponential tendon stiffness model (see Eq (10)). The exponential tendon model describes a linear relationship between stiffness and force within the toe region as observed by Rack and Westbury [2]. The polynomial model as described by Buchanan does not have this property. Therefore, the exponential model is used throughout this work.

Tendon model with exponential force length relation by Rack and Westbury.

A linear force stiffness relationship, as described by Rack and Westbury [2], can be reached by an exponential force strain model. In general the stiffness of a material is defined as(4)

Rack and Westbury could show, based on measurements of the cat soleus tendon, that the stiffness of a mammalian tendon k_T depends linearly on the applied force [2].(5)

Proske and Morgan concluded that the linear stiffness behaviour of Eq (5) is only valid in the toe region and assumed that the tendon stiffness beyond the toe region is constant [20], thus leading to a linear force strain behaviour beyond the toe region (Fig 2A). Based on these findings, the following differential equation models the stiffness for during the lengthening of the tendon (Δl_T ≥ 0).(6)

Solving the differential equation for each case leads to the following expression, if c₂ ≠ 0 is assumed.(7)Where Δl_toe is the change in tendon length at which the toe region ends (see Fig 2A). A₁ and A₂ are coordinates in the solution space of Eq (7), which also contains physically implausible solutions. In order to gain a physically feasible solution, the condition from Eq (2) is applied and continuity of for Δl_T > 0 is assumed. This results in the following, unique solution.(8) (9)

The substitution of A₁ and A₂ in Eq (7) leads to the following model description.(10)

Here, the tendon parameters are included explicitly in the function definition as parameter vector .(11) Eq (10) is not defined for compression (Δl_T < 0). It is widely assumed that the tendon resists with nearly no force to compression [19, 21]. This compression behaviour, however, is not needed in this study as the tendon will only be observed in extension. Due to the type of data acquired in this study, the inverse of Eq (10) is needed.(12) Eq (12) is used for the tendon model that was used to represent the experimental data as described in the following section.

Methods for tendon modelling based on in vivo measurements

Experiments were conducted with 5 subjects from which 3 are shown in this study. All subjects where recruited by call for participation and gave written informed consent to the pseudo anonymous usage of the recorded data. All subjects that responded to the call for participation and stated that they had no related injuries where accepted for the study. The study was conducted according to the ethical guidelines of the German Society for Psychology (DGPs) and the German Psychologists Association (BdP), and approved by the Ethics Committee of the University of Bielefeld (EUB 2017–156 02.08.2017). The experiment data is published in a zenodo data repository [22].

The displacement of the lower biceps brachii myotendinous junction ΔL_mt was measured over 7 elbow angles Θ and over 6 levels of desired force F_w for each subject. The force was measured at the subjects wrist in isometric conditions.(13)

During the experiment, the right arm of a subject was positioned on a support rack, as shown in Fig 3. The support rack was used to keep the elbow angle constant during the experiment. A force sensor consisting of a load cell and strain gauges in half bridge configuration was used to measure the force at the wrist , while the hand was fully supinated. A hook-and-loop strap mechanically connected the wrist to the load cell. The lengthening of the wrist strap was assumed to be negligible in terms of changes in Θ (linear stiffness of the wrist strap is at a preload force of 15N). The strain gauges were processed with a bridge measurement module (NI-9237, National Instruments Corp., Austin, TX, USA) and was logged via LabVIEW (National Instruments Corp.) at a sampling frequency of f_F = 1kHz. The measured force was offset by the wrist strap preload at the beginning of the experiment, filtered with a moving average filter over 40 samples and fed back to the subject via a virtual gauge displayed on a computer screen using LabVIEW.

[Figure omitted. See PDF.]

The distance l_w(Θ) between the elbow and the wrist strap, which mechanically connects wrist and load cell, is marked with gray lines. The distance l_w depends on the elbow angle Θ, which is shown in black. A change in joint angle is constrained by the depicted construction of aluminium profiles. If the screw in the hinge joint of the profile structure is loosened, the joint angle can be adjusted.

The myotendinous junction was observed using an ultrasonic device in B-mode. The ultrasonic transducer was handled by an experienced operator. As indicated in Fig 3, the tendon displacement together with the wrist force was measured in three steps.

1. In the unloaded muscle state (no wrist force) a video capture of the ultrasonic B-mode image was started. (t₀ = 0s)

2. The subject was instructed to flex the right arm such that the desired force level could be measured at the wrist (feedback via computer screen). (t₁ = 2s)

3. The B-mode video capture was stopped and then the subject was instructed to relax again. (t₂ = 8.5s)

4. At least 20 seconds of rest.

These steps were initiated by beeping sounds at the given times t and repeated for each desired force level F_w at each elbow angle Θ. For each subject the experiment started at an elbow angle of Θ = 90^∘. First, the maximal voluntary contraction force was measured without ultrasonography imaging. Then each desired force level was measured three times in ascending order before the elbow angle was reduced by 15^∘ (note Eq (13)).

Later, the movement of the myotendinous junction was extracted from the B-mode video captures.

Additionally, distances between anatomic landmarks of the arm were measured with a tailor tape measure (S1 Table). These distances were used as the basis of a skeletal arm model (cmp. Eq (37)). With this model, the tendon properties were determined by means of a fitting approach. Here, it was assumed, that any observed movement of the myotendinous junction parallel to the humerus is attributed to the lengthening of the lower tendon of the biceps brachii muscle during contraction. Using force and shortening information, a local tendon stiffness for a given elbow angle Θ can be estimated or a global characteristic curve which is valid for several elbow angles can be constructed. The latter was the goal of this work.

In vivo measurement of the tendon displacement

In order to quantify the displacement of the lower biceps brachii myotendinous junction, ultrasonic video images of the anatomical structure were recorded with an E1 ultrasound device equipped with the 50 mm linear ultrasonic array transducer L741 (SonoScape Medical Corp. Shen Zhen, P.R. China) in B-mode at a frequency spectrum of 9.5MHz to 12.2MHz. The linear array was chosen as its field of view allows safe observations of the myotendinous junction during the complete shortening process up to the maximal voluntary contraction. The intensity in a B-mode image corresponds to the amplitude of reflected acoustic waves. Reflections occur at tissue changes due to varying material properties. A homogeneous tissue section would appear in black, in contrary to inhomogeneous tissue sections or junctions between two homogeneous sections. Fig 4A shows the lateral view of the measurement setup as an anatomical sketch with an orientation of the linear array (SP) that corresponds to the imaging orientation in Fig 4B (left, marked I). In Fig 4B the actual two ultrasonic images of the myotendinous junction in two imaging orientations I (Fig 4B, left) and II (Fig 4B, middle) are given. To further illustrate the position and orientation (imaging plane I or II) of the linear ultrasonic array, an anatomical sketch of the biceps brachii and brachialis from an anterior perspective is given on the right side of Fig 4B. Imaging plane I (left side of Fig 4B) represents the parasagittal plane crossing the humerus (Hum) and the aponeurosis (BicAP). The aponeurosis connects the tendon (BicT) to the muscle (cmp. Fig 4A). Here, each muscle fiber’s terminal Z disk is connected to the tendon collagen structure [23]. Due to the change in material properties, the aponeurosis appears as a distinct feature in the ultrasonic image. As the muscle fibers gradually terminate, the muscle belly flattens out. A distinct region is identified where the aponeurosis emerges from the main tendon structure (furcation of BicT ranging into the biceps Bic). This point is referred to as , and is marked by a white cross in Figs 4B and 5. Contraction of the muscle leads to movement of that point. Due to the anisotropic structure of skeletal muscle, the biceps brachii can be discriminated from the brachialis muscle group (Bra). The cross-sectional view in imaging plane II allows a clear distinction of the individual muscles, as the the muscle structure in the abaxial plane is mostly homogeneous as opposed to the tissue transitions.

[Figure omitted. See PDF.]

In (A) the lateral view of the experimental setup in form of an anatomical sketch is shown. An ultrasonic probe (SP) is used to image a cross-section of the distal part of the upper arm. The cross-section contains the distal tendon of the biceps brachii (BicT). Simultaneously, the force generated by the arm is measured with a force sensor (FS) near the wrist (cmp. Fig 3). The biceps brachii (Bic), its long head (BicLH), its short head (BicSH), brachialis (Bra) and humerus (Hum) are also labeled. In (B) the ultrasonic images of two orthogonal cross-sections are depicted. The orientation of both imaging planes are marked in a sketch on the right side, which shows the anterior perspective on the Bic and Bra. Orientation I was the only imaging plane used during the presented study (cmp. (A)). The structures of the Bic are marked in red. The junction of Bic and BicT was observable in this view (white cross) as well as part of the aponeurosis (BicAP). Imaging plane II is orthogonal to imaging plane I. The furcation (wing-like structure) of the aponeurosis was observable (marked in orange) as well as different parts of the Bra muscle group—the brachialis deep head (BraDH) in green and the brachialis superficial head (BraSH) in blue.

[Figure omitted. See PDF.]

(A) shows the unloaded state n_u (F_w = 0) and (B) the loaded state n_l. Reference points on the humerus and on the myotendinous junction are marked in the frame of each state. The reference points were used to calculate the position of the myotendinous junction relative to the bone local coordinate system (brown). To calculate the bone surface orientation within the frame, a second marker was placed arbitrarily on the bone surface. That is for the unloaded state and for the loaded state. The bone line was set such that it included both bone points and . The bone line is drawn in turquoise for each frame.

During the experiment the ultrasonography operator focused on maintaining the imaging plane I (cmp. Fig 4B). Hence, the operator had to constantly compensate changes of the probe orientation, induced by skin surface changes during contraction. The operator focused on maintaining the same reference points marked by brachialis insertion into the humerus, as well as the myotendionous junction within the frame. If this was not possible, the operator instructed a repetition of the measurement.

Manual extraction of the myotendinous junction displacement.

In Fig 5 two B-mode frames in the parasagittal plane are shown (for reference compare with orientation of imaging plane I in Fig 4A). These frames were selected by an experienced ultrasonography operator who decided, based on the movement of the anatomical structures, when the steady states of the myotendinous junction movements were reached. The first steady state was the unloaded state whereas the second state was reached when the subject applied the instructed force level. Due to deformation of the tissue during contraction, the two frames show slightly different viewpoints. During the experiment, the ultrasonography operator chose a fixed reference point and made sure that it stayed within the frame. As reference points, insertion points of brachialis fibres at the humerus were chosen as these are fixed within the coordinate system of the humerus. These points served as reference coordinate system. Fig 5 depicts the bone coordinate frame at that reference point in brown for the unloaded Fig 5A and loaded state Fig 5B. It is oriented such that the abscissa is congruent with the bone surface. The bone surface in the B-mode image is abstracted as a line (Fig 5, drawn in turquoise). The bone line was parametrized by means of of the orientation angle (see Fig 6) and the bone reference point (where the tilde symbol indicates that the respective quantity is linked to the manual data extraction). The orientation angle was derived from the bone reference point and a second marker . These parameters were manually marked for the unloaded state n_u and loaded state n_l.

[Figure omitted. See PDF.]

A sketch which demonstrates how the bone line angle was calculated during manual data extraction () and automatic data extraction (). (A) demonstrates the case of manual bone line labeling in which was derived from the vector as described in Eq (16). Analog (B) demonstrates how was derived from a bone line expressed in Hough-line parameters (cmp. Fig 9) as formulated in Eq (24).

The point of the myotendinous junction was chosen where the aponeurosis crosses with an imaginary extension of the fascia in both states. The junction point was transformed into the reference coordinate system for each frame.(14)Where is the bone line orientation angle relative to the abscissa of the bone local coordinate system. And R(•) is the rotation matrix for the bone local coordinate system. The bone local coordinate system is a left handed coordinate system. rotates in the clock wise direction for positive angles. Consequently R⁻¹(•) and are described as follows.(15) (16)Where is the vector pointing from the bone reference point to the bone orientation marker as shown in Fig 6A. The operator [•]_x returns the x coordinate and [•]_y the y coordinate of a vector.

Calculation of myotendinous junction movement.

After the myotendinous junction was described in the bone local coordinate system, the tendon lengthening could be calculated for each frame. If the absolute position of the bone reference point in a skeleton model and the path of the tendon were known, the lengthening of the tendon could be calculated accurately. In this work, the lengthening of the tendon was assumed as the movement of the myotendinous junction from unloaded to the loaded state parallel to the bone line.(17) (18)

Data for of subject 0 is presented as an example in Fig 7A. Later the quality of the extracted data is discussed and compared with a (semi-) automatic approach.

[Figure omitted. See PDF.]

The length of the junction movement parallel to the bone of subject 0 is shown for the manually extracted data as depicted in (A) and for the automatically extracted data in (B). In both cases the mytendinous junction movement parallel to the bone is shown over the wrist force F_w for all elbow angles Θ (cmp. Eq (13)). The bone line angle as well as the evaluated frames are adopted from the manual evaluation as described in Eq (31). In (C) the difference between the automatically extracted data and the manually extracted data of subject 0 is shown.

Computer aided extraction of the myotendinous junction displacement.

In addition to the manual procedure for extracting the myotendinous displacement (see Eqs (17) and (18)) a computer-aided procedure was employed. The computer-assisted method allows a) to minimize any inaccuracies due to the experimenter’s subjective judgment during manual extraction of reference points, b) to ensure reproducibility of results, and c) to achieve higher throughput when processing large data sets.

In order to calculate the junction movement , the computer-aided method extracts the myotendinous junction , the relative orientation of the bone and the bone reference point from B-mode ultrasound frames based on different algorithmic methods as depicted in Fig 8 and described in the following paragraphs. The hat symbol is used as annotation of measurement quantities related to the computer-aided feature extraction.

[Figure omitted. See PDF.]

The computed features are from left to right: the mytendinous junction, bone line and bone reference point. The mytendinous junction is tracked with a neural network based on the MobileNetV2 algorithm. The bone line is tracked via Hough-Line-Transformation. The bone reference point is tracked by a method based on the Lucas-Kanade algorithm.

Detection of myotendinous junction using deep neural networks.

In order to detect the myotendinous junction from the data set of B-mode ultrasound frames, a deep learning approach using transfer learning was used (see Fig 8). Automated tracking of the myotendinous junction using convolutional neural networks has recently been shown to provide accurate results comparable to manual tracking methods [24, 25] and to be less error-prone than semi- and fully-automated myotendinous junction tracking approaches based on optic flow [26, 27] or block-matching algorithms [28, 29].

To generate and train a convolutional neural network capable of detecting the myotendinous junction of the biceps brachii from ultrasound images, the markerless pose-estimation framework DeepLabCut [30] was used. The architecture of the underlying convolutional network is based on MobileNetV2 [31] and was pretrained on the ImageNet data set by [32]. A training data set was generated via random selection of 203 video sequences from the data set of B-mode ultrasound images. From each individual video sequence 10 images were extracted in a randomly and temporally uniformly distributed manner resulting in a total of (n =) 2030 frames. In each of the extracted image frames, the myotendinous junction of the biceps brachii was labeled by hand, hence, serving as a ground truth. The labeled images were subsequently split into two data sets for use in training (n_t = 2185 frames) and evaluation (n_e = 115 frames). Training of the neural network was performed using the training data set with a batch size of 20 for 200000 iterations and DeepLabCut’s standard parameters. The performance of the neural network was subsequently evaluated via computation of the the mean average error (MAE) of the Euclidean distance between the predicted and the ground truth location of the myotendinous junction in the evaluation data set. This resulted in an MAE of 4.32 pixels corresponding to 0.50mm. After training, the neural network was used to automatically infer the location of the myotendinous junction from the individual video sequences of B-mode ultrasound images (hat symbol indicates computer-aided feature extraction).

Computation of the bone orientation.

The bone orientation relative to the x-axis was previously calculated from two manually selected points on the bone surface . Due to the lack of explicitly definable features on the bone surface for these points in the B-mode frames, feature tracking analog to the myotendinous junction was not an option. The bone surface appears as a continuous feature in the B-mode image as depicted in Fig 5. However, the bone surface has a high contrast as compared to the surrounding muscle tissue. The transition from muscle tissue to bone can be abstracted as a line. Therefore, the B-mode image was filtered for edges and a line detection algorithm was applied. As the fibers in the muscle tissue also have a line characteristic, the B-mode image was cropped such that the risk of a line feature in the muscle tissue being detected as bone line was minimized. A block diagram of this process is shown in the middle section of Fig 8. In order to confine the algorithm to an image region relevant for the detection of the visual line separating muscle from bone tissue, each ultrasound image frame was cropped using the FFMPEG library [33] resulting in a size of (w_c =) 600 by (h_c =) 414 pixels of the cropped frame. (x_c = 320, y_c = 350) is the horizontal and vertical pixel position of the top left corner of the cropped region within the original frame (see Fig 9A). After cropping, the Canny edge detection algorithm ([34], as implemented in OpenCV [35]), was applied to each frame. This algorithm contains a four-stage process in which 2 of the 4 stages need to be parameterized. First, a noise reduction with a 5x5 Gaussian filter is applied. Second, the intensity gradient is computed. This process is parameterized by the aperture size (a_s = 3) and the norm used in gradient computation. Here, the computationally efficient L₁ norm was used. Third, after gradient computation, a non-maximum suppression is applied. The result is an image in which the value of each pixel describes the gradient of an edge. Finally, in the fourth step, these edges are filtered by hysteresis thresholding. In this work, pixel with a gradient value above grad_max = 220 were accepted as edges. Pixel with a gradient value above grad_min = 50 were considered edges if they were connected to edges above grad_max. All other pixel/edges were filtered out. The used parameters are shown in S2 Table.

[Figure omitted. See PDF.]

(A) is a frame of the cropped section of a B-mode video. The detected lines have the same color coding as in the scatter plot of the bone line parameter space. (B) The corresponding scatter plot of detected Hough-Lines in their parameter space. The biggest detected cluster is shown in yellow and its mean is shown in red. All other data points are black.

The remaining edges were used during the Hough-transformation for line detection. The Hough-lines are defined by a vector in cylindrical coordinates (r, φ). The vector has to be understood as a point on the described line . Where is the parameter form of the bone line, which is defined as follows.(19)

Furthermore, it is assumed that is oriented orthogonal to that line (). Because of these two definitions, the bone line l_b(h) is fully defined by the tuple (r, φ) where r is in pixel and φ in radiant. These parameters span a constrained parameter space of possible lines.(20)

This Hough-line parameter space was transformed into a constrained and discretized Hough-line search space which is described by(21)and(22)where Δr = 1 and are the search resolution parameters.

For each edge pixel in a frame all lines containing this pixel in the discretized parameter space were considered and the counter of the corresponding parameter tuple was increased. After all edge points had been processed, those parameter tuples with counter values larger than a threshold (Threshold_Hough = 40) were considered to be bone line candidates. These candidates were sorted by their counter values. The candidates with the 10 highest counter values were selected for further processing (see Fig 9A). The parameters used for the OpenCV HougLines function were shown in S3 Table in the S3 Table. If a function parameter is not listed, the default value was chosen. The selected Hough-lines were used to calculate a mean bone line. As some structures in the muscle tissue also resulted in detected lines, a clustering algorithm (DBSCAN from SciPy [36]) was applied to filter for bone lines, as shown in Fig 9B. By means of majority voting, the biggest cluster was chosen and its mean parameter vector was used to describe the bone line. In subsequent steps, only the orientation of the bone line was used. Before clustering, the Hough-line parameters were normalized on the expected range for a bone line.(23)

This way distances between parameters in and can be treated equally, which is beneficial for clustering. The function parameter m ∈ {0, 1, …, 9} is the index of a bone line parameter set in frame n. After this normalization, DBSCAN is applied which is a density based clustering algorithm defined by [37]. There are 3 types of data points. 1. a core point, which has at least n_min,s − 1 samples within the distance ϵ_core. 2. a border point, which is within the distance ϵ_core of a core point. 3. noise points / outlier. The latter are not reachable from any point. A point is called reachable, if there is a chain of core points within the distance ϵ_core. Points, which reach each other, are part of the same cluster. For bone line clustering, the amount of core points was set to n_min,s = 2 and ϵ_core = 0.02 was chosen. As distance function, the Euclidean-distance L₂(•) = ‖•‖₂ was used. In S4 Table the parameters are listed with their names as given in the SciPy module [36]. Default values were used for those function parameters that are not listed. The bone line was then calculated as the mean of the biggest cluster (r_b(n), φ_b(n)) where n is the index of the corresponding frame. Based on this findings, the bone line orientation angle relative to the abscissa was calculated as follows.(24)

A derivation of this relationship is shown in Fig 6B. In the used coordinate system, a clock wise rotation is positive. This description of can be used instead of the manual labeled bone orientation .

Tracking of the bone reference point via Lucas-Kanade optic flow estimation.

As compared to the myotendinous junction, the visual texture of bones in B-mode ultrasound images does not provide distinct visual features consistent for all subjects (Fig 4). Hence, the automated inference of reference points using a deep neural network approach based on transfer learning—such as described in the section Detection of myotendinous junction using deep neural networks—is not feasible, as such an approach would require supervised training based on individual features of the bone surface for corresponding subsets of ultrasound images. Instead, for the automated tracking of the bone reference point on the humerus from sequences of ultrasound images, an optic flow based approach was used. Lucas-Kanade optic flow estimation [38] is based on the assumption that the distribution of apparent velocities induced by motion is constant in a local neighbourhood of an image pixel under consideration. This allows for the temporal tracking of visual features in a sequence of images which do not share distinct textural properties across a set of image sequences. Recent studies have used the Lucas-Kanade based estimation of optic flow to temporally track visual features in ultrasound images [26, 27].

To obtain the bone reference point from the data set of B-mode ultrasound images, a region of interest (ROI) representing the bone structure of the humerus was computed for each frame of the data set of ultrasound images. The ROI was defined as the area below the visual line separating muscle from bone tissue within the cropped image frames (see paragraph above; Fig 9A). Within the ROI, visual features were selected for the first frame of each ultrasound image sequence using Shi-Tomasi corner detection [39]; implemented using OpenCV [35]. The maximum number of corners returned via the corner detector was set to N_c,max = 1000 while corners with a quality level of Q_c < 0.25 and a Euclidean distance of d_c,max < 3 image pixels were rejected. The size of an average block for computing a derivative covariance matrix over each pixel neighborhood was set to s_b = 3. This resulted in a mean of corners (standard deviation σ_c = 24.48) detected per analyzed image frame. Subsequently, the detected visual features (or corners, respectively) in the first frame of an image sequence were iteratively tracked in the remaining frames of the corresponding image sequence using a pyramidal implementation of the Lucas-Kanade feature tracking algorithm [40]; implemented using OpenCV [35]. The maximum number of pyramidal levels taken into account was set to N_p,max = 4 and the size of the search window at each pyramid level to . The termination criteria of the iterative search algorithm were set to ϵ_p = 0.03 for the desired change in parameters at which the algorithm stops and i_p,max = 30 for the maximum number of iterations to compute. The bone reference point was subsequently obtained for each frame n in each sequence of ultrasound images via computation of the centroid (i.e. the mean) of the visual features with coordinates [f_x, f_y]^T tracked via Lucas-Kanade optic flow estimation according to(25)whereas features leaving the ROI during tracking within an image sequence due to movement of the bone in relation to the ultrasound sensor were discarded.

Comparison of manual and automatic feature extraction.

The automatic feature extraction was compared to the manually extracted data per feature.

When evaluating the tracking of the myotendinous junction, the absolute movement between states (unloaded, loaded) was compared,(26)as the investigator in the manual method and the automatic method did not necessarily track the same feature of the myotendinous junction. Tracking of the bone line was evaluated separately for the loaded and unloaded state as the angle of the bone line was comparable on a frame basis.(27)

In case of the bone reference point the absolute movement of a manually chosen feature on the bone (see Fig 5) and the absolute movement of an automatically chosen feature by Lucas-Kanade was compared. The comparison of the absolute distance between steady states(28)

The results of these comparisons are shown in Fig 11. For each comparison a mean absolute error (MAE) and the Pearson correlation coefficient (R) is calculated.

Calculation of myotendinous junction displacement.

The myotendinous junction displacement parallel to the bone needed to be calculated analog to the manual data extraction. For that, Eqs (14) and (17) still held true. For each manually extracted feature (, , ) there is an automatically extracted counterpart (, , ). These features are compared in Fig 11. It appears that the error between manually and automatically labeled absolute bone reference point movement (see Eq (8)) correlates with the wrist force F_w (R = −0.58). This finding is later discussed. Therefore, for further steps the bone reference point was used from the manually labeled data set instead. Consequently, the myotendinous junction displacement vector is defined as follows.(29) (30)

Then the displacement parallel to the bone is defined as:(31)

For the steady states n_u and n_l, the same frames were picked as for the manual data extraction. Data for is shown in Fig 7B for subject 0, as example.

Measurement of elbow joint torque

For the prediction of the myotendinous junction movement, the tendon force needs to be calculated. During the experiment, the wrist force F_w was measured. From this force, the net elbow torque was calculated for the unloaded and loaded state. Later in this work, this data was used to estimate the torque which is generated by the biceps brachii from which the tendon force can be estimated.

Wrist force extraction.

The unfiltered wrist force and data from B-mode ultrasound frames were not synchronized. However, steady states can be identified in . Since it was unknown at which time during the steady state the frames were captured, the mean of the wrist force in each steady state was calculated. The force contributing to a moment at the elbow joint during the unloaded state was assumed to be zero. A measured force during this time window therefore had to be attributed to the preloading of the wrist strap. The measured mean force during the loaded state was therefore offset by the mean force during the unloaded state in order to obtain the actual elbow moment generating force F_w.(32)Where n_l0 is the first sample of the loaded state and n_l1 is the last sample of the loaded state. The same principle applies to the unloaded state (n_u0 first sample—n_u1 last sample).

Torque at the elbow joint.

The moment arm l_w(Θ) changes over the elbow angle due to the construction of the support rack (Fig 3). For this reason, the distance between the elbow joint and the wrist strap was measured for each elbow angle Θ. With this measurement, the corresponding elbow joint torque τ_e was calculated.(33)

This elbow torque was later used to estimate the torque which was generated by the biceps brachii.

Modeling the tendon and parameter fitting

As described in the section Mammalian tendon modeling, there are two major tendon models. In this work, the model was used. It is derived from elastic properties as found in experiments measuring the stiffness of a mammalian tendon in vitro. As the experiment in this work was conducted in vivo, the variables l_0,T,l_T, , F_T were not known a priory.

The absolute tendon length l_T and the resting tendon length l_0,T remained unknown, whereas the tendon force F_T and the stiffness parameter vector have been estimated. The tendon force was estimated by means of a two dimensional tendon path model. The stiffness parameter vector was estimated by fitting the tendon model on the measurements of the myotendinous junction displacement. Instead of the absolute tendon length, the change in length from rest to the current state was used Δl_T = l_T − l_0,T. It was estimated based on the myotendinous displacement parallel to the bone (as observed in the B-mode videos) and by the assumption:(34)

In this study, only the distal biceps brachii tendon was observed with ultrasonography. The force exerted on the distal tendon was estimated based on the measured elbow joint torque τ_e as given in Eq (33) and a model of the abstracted tendon path (cmp. Fig 10A). The tendon path model was parameterized by measurements of anatomic landmarks as given in Eq (37).

[Figure omitted. See PDF.]

(A) Schematic depiction of the applied biomechanical model of muscle and tendon in the elbow joint. k_T1 and k_T2 are non linear spring elements representing the lower and upper biceps brachii tendons. The shoulder and lower arm insertions are also marked ( and ). These points were determined based on anatomic landmarks as formulated in Eq (37). The elbow joint was abstracted as revolute joint. It was located at . The effective lever length of subject 0 is shown in (B) over the elbow joint angle Θ. The maximal effective lever length is marked at (Θ = 93.3^∘, l_eff = 5.2cm).

The parameter vector was estimated with a fit on the collected subjects’ data. As a target function for the fit, the myotendinous displacement ΔL_mt(F_w) as a function of the wrist force was used. This is because the wrist forces F_w were voluntarily adjusted by the subjects during the experiments and have only a small variance per force interval (see Fig 7). For the myotendinous displacement ΔL_mt(F_w), the formulation for from Eq (12) was used.

Estimation of the tendon force.

In this work, tendon displacement was predicted based on the measured wrist force F_w and an abstraction of the skeletal geometry as shown in Fig 10A. The model results in an effective lever l_eff. From this, the tendon force could be calculated by(35)since the elbow joint was assumed to be a revolute joint. The torque τ_bic is generated by the biceps brachii. It was estimated based on the elbow joint torque, as it was not directly measured.(36)

In this equation, τ_r is the torque generated from other sources (i.e. friction, other muscle groups, surrounding tissue). The forearm, which contains two bones (ulna and radius) as mechanical base structure, was abstracted as a single lever. The distal biceps brachii tendon T₁ inserts into this lever at with the same distance l_i. As distance l_i the distance between the elbow and the insertion into the radius was used. The fascial insertion of the distal biceps brachii aponeurosis was ignored.

The two headed biceps brachii was abstracted to a single head inserting at into the shoulder. With the coordinate systems origin placed within the elbow revolute joint, the points(37)were defined (see Fig 10A). Here, the distance from the distal tendon insertion point to the elbow joint l_i was assumed to be proportional to the length of the ulna (l_i ∝ l_u).(38)

The conversion factor was set based on the proportions depicted in figure. 334 of [41]. The effective lever length can be calculated as follows.(39)Where is a unit vector in the direction of . The course of the effective lever length l_eff over the elbow angle for subject 0 is shown exemplarily in Fig 10B.

With Eqs (35) and (39) the estimation of tendon force is complete. It was used for fitting the tendon model onto the extracted myotendinous junction movement and the wrist force F_w.

Fitting the tendon model to data.

Due to the small amount of data per elbow angle and its high variance in the myotendinous displacement ΔL_mt, it was preferable to fit the model onto all angles in one process instead of breaking it down into several sub-processes—one for each angle. The material properties of the tendon were assumed invariant, also in terms of changes of the elbow angle. However, if the predicted biceps torque is less than the elbow torque (τ_r > 0, see Eq (36)) the tendon parameters are effected. This means that, when τ_r depends on the elbow angle, it seems as if the tendon parameters are changing depending on the angle. The same holds true for the effective moment arm l_eff. These effects are further analyzed in the results and evaluation section. Therefore, it was assumed that only a subset (or sub-range) of the measured elbow angles can be explained by the current model. Within this sub-range, τ_r ≈ 0 was assumed. Additionally, l_eff was assumed to predict the effective lever length for that sub-range. As the valid elbow angle range was not known a priori, it was narrowed down in an optimization process. In this optimization process, fits of the model for decreasing ranges of Θ were performed. Later, the elbow angle range was used as a fit result whose mean cost per sub-range is stable as compared to all angle ranges (cmp. Figs 12A and 13A, example data from three subjects). For each such elbow angle range fit, the error function(40)was minimized. The error function g calculates a sum of differences between the predicted and measured tendon displacement. The function μ(•) was used to distort the error function to cope with variance in the data. In this case, μ(•) = ln(•+ 1) was chosen. The functions Θ[k], F_w[k] and are discrete functions returning the respective values of a data point indexed by k. The set (lb for lower bound) contains all indices of data points where Θ_lb ≤ Θ. Therefore, it is possible to calculate fits over subsets of elbow angles.(41)

This minimization was done by performing fits with(42)as lower boundaries for the elbow angle subsets which resulted in parameter sets for each value of Θ_lb. The mean cost per subset was calculated as follows.(43)Where is the amount of data-points of the corresponding angle subset defined by Θ_lb. As an optimization algorithm, the non linear function optimizer trust region reflective from the Python SciPy package [36] was used. For the optimization process(44)was used as the initial parameter vector. As lower and upper boundaries(45)were used. Here, eps₆₄ is the machine epsilon of the 64 bit float datatype. Results of the fitting process are depicted in Figs 12 and 13 and evaluated in the following section.

Results

In this work, a manual and a (semi-) automatic method of extracting data (myotendinous junction displacement in a bone local coordinate system) from B-mode ultrasound images were used. Both methods where evaluated based on the displacement of the distal biceps brachii tendon during the generation of different wrist forces at different elbow angles. Subsequently, tendon models were fitted to the manual and (semi-) automatically extracted data from three subjects (cmp. Figs 12 and 13). Special emphasis was placed on fitting the models not only to one elbow angle but to subsets from the entire range of elbow angles.

Comparison of manual and (semi-) automatically extracted data from B-mode ultrasound videos

The myotendinous displacement parallel to the bone was either extracted manually or automatically from frames of B-mode ultrasound videos. The two methods are compared for subject 0 in Figs 7 and 11. During automatic feature extraction, for each feature (myotendinous junction, bone line and bone reference point) a distinct algorithm was applied (Fig 8: left, middle and right, respectively). Thus, each of the three algorithms needed to be evaluated separately. For each feature tracking the result was compared with the manual tracking as ground truth. The difference from manual to automatic evaluation is shown over the wrist force in Figs 7C and 11.

[Figure omitted. See PDF.]

The green check mark in the upper right corner indicates that the corresponding method was used for evaluation in this work. The red cross indicates that the manually labeled data was preferred over the automatically extracted data.

Tracking of the myotendinous junction.

The comparison of the absolute myotendinous junction displacement between manual and automatic data extraction resulted in a mean absolute error (MAE) of 0.67 mm ≐ 9:0 pixel. The differences were distributed approximately symmetrical around zero, with two outliers. As for the outliers, the neural network extracted a larger movement than the manual evaluation. The difference of absolute movements correlates with the wrist force by R = −0.07 (see Fig 11, left), where R is the Pearson correlation coefficient.

Tracking of the bone line.

Tracking of the bone line was evaluated separately for the loaded and unloaded state. In both states, the computation of the mean absolute error resulted in MAE < 1^∘ (see Fig 11, middle), which was below the resolution of the Hough-line transformation. Apart from outliers, the mean absolute errors for both states were distributed around zero. This is due to image cropping, which mostly affected the unloaded state due to a smaller distance between bone and ultrasonic probe. Using a higher cropping window would have led to false bone line detection within muscle tissue. In either, the loaded and the unloaded state, the bone line angle γ correlated with the wrist force F_w by R = −0.08 (Pearson correlation coefficient).

Tracking of the bone reference point.

The comparison of the absolute distance between steady states correlated with F_w by R = −0.58 (Pearson correlation coefficient). This resulted in an error distribution around an offset value unequal to zero.

Comparing the (semi-) automatic data extraction with the manual data extraction for subject 0 led to a MAE of 1.15mm (see Fig 7C). Therefore, with a mean absolute error difference close to one millimeter, the results of both methods were assumed to be equivalent.

Results of the tendon model fitting process

As described in the Modeling the tendon and parameter fitting section, elbow angle sub-ranges of the manually and (semi-) automatically extracted data were used to fit the tendon model for three subjects as examples. The fitting results of the manual and (semi-) automatically extracted data are shown in Figs 12 and 13, respectively.

[Figure omitted. See PDF.]

(A) Myotendinous junction displacement as manually extracted from B-Mode videos plotted over the wrist force F_w. (B) Normalized costs of each angular subset fit (see Eq (43)). (C) Fits of the angular subsetsΘ_lb = 30^∘, Θ_lb = 45^∘,Θ_lb = 60^∘,Θ_lb = 75^∘. (Same color coding as in Fig 7).

[Figure omitted. See PDF.]

(A) Myotendinous junction displacement which is (semi-) automatically extracted from B-Mode videos plotted over the wrist force F_w. (B) Normalized cost of each angular subset fit (see Eq (43)). (C) Fits of the angular subsetsΘ_lb = 30^∘, Θ_lb = 45^∘,Θ_lb = 60^∘,Θ_lb = 75^∘. (Same color coding as in Fig 7).

In case of the manually extracted data, the fits over elbow angle sub-ranges showed that there are sub-ranges with nearly identical mean cost ϵ (cmp. Fig 12B). For smaller Θ_lb (larger angle sub-range), the mean costs increased for all three subjects.

In case of automatically extracted data, this trend only held for subject 0 (cmp. Fig 13B).

For subject 0 and 1 the optimal elbow angle sub-range was identical for (manual) and (automatic). For subject 2 the optimal subset deviated between manual (Θ_lb = 75) and (semi-) automatically extracted data (Θ_lb = 30).

It was common for all fits that with decreasing angle, a higher tendon displacement was predicted. In contrast, the measured tendon displacement generally decreased for each elbow torque if the elbow angle became smaller. Thus, the model had a tendency to overestimate the myotendinous junction movement for decreasing elbow angles.

In order to compare these fitting results with the findings in [9] a linear stiffness was calculated for each fit on the angle subranges Θ ≥ 30 (see Table 1). The linear stiffness was computed analog to [9] with the tendon lengthening at 20%Â max(F_T1) to 80% max(F_T1). An example linear stiffness is shown in Fig 14B in comparison to the nonlinear tendon stiffness function.

[Figure omitted. See PDF.]

(A) Validation of the maximal effective lever length. The maximal effective lever length l_eff over the radius length is shown as black dots, as presented in the cadaver study by Murray et al. [42]. The linear fit of Murray’s data is shown as black line. The blue stars are the maximal effective lever lengths of the subjects in this study. A complete course of the effective lever length is shown in Fig 10B for subject 0. (B) Non linear and linear prediction of myotendinous junction movement. The non linear prediction of the myotendinous junction movement over tendon force is shown in orange. The prediction is based on Eq (10). As parameter set, the resulting fit on (semi-) automatically extracted data from all angels Θ ≥ 75^∘ of subject 0 was used. In addition, a linear model is shown in black. This linear model was parameterized such that the linear function includes the predicted values for 20% and 80% of the maximal tendon force generated during voluntary contraction . The stiffness of the depicted linear model was c_lin,20-80 = 362N/mm.

[Figure omitted. See PDF.]

In addition, the optimal angle sub-range is highlighted in gray for manually and automatically extracted tendon displacement (cmp. Figs 12B and 13B).

Discussion

In the work described here the displacement of the biceps brachii myotendinous junction parallel to the humerus was extracted from B-mode ultrasound images as part of an experimental setup that also allowed determination of applied loading forces and elbow joint angles. A manual and an automatic method to extract the relative motion of the myotendinous junction based on visual features from ultrasound images were developed. The performance of both methods was compared. Using experimental data from three subjects, a tendon model was fitted to the data of both methods used for extracting the myotendinous junction displacement. In both cases a reliable fit was achieved.

Experimental setup

In the experiment, the myotendinous junction was observed with a linear ultrasonic probe array, while the wrist force was measured asynchronously with respect to the ultrasonic measurement (cmp. Fig 3). As a consequence, the collected data and the resulting fits could only be used for the examination of steady states of the tendon behaviour. Any dynamic processes such as creepage were not taken into account.

The myotendinous junction was observed in a saggital plane with a linear ultrasonic probe array (cmp. view I in Fig 4B). The moytendinous junction, and also the course of the aponeurosis, has a three-dimensional shape and was not oriented in a particularly aligned way in the saggital plane (cmp. Fig 4B). Therefore, during muscle contraction, the motion of the myotendinous junction had components within the saggital imaging plane as well as outside the plane, e.g. perpendicular to it. The latter were not detected in the presented setup.

Additionally, the biceps muscle bulges during contraction, changing the outer contour of the arm, which also changed the relative orientation of the sensor probe to the myotendinous junction. As the ultrasonic sensor probe was handheld, the operator attempted to compensate tilting during the execution of the experiment. If the sensor probe tilted too much the bone reference point moved out of the imaging plane and therefore could not be tracked. These cases were excluded from the data set.

A future extension of the setup could include a two-dimensional ultrasonic sensor array, which would support the investigation of a three-dimensional measurement volume. This would allow consistent tracking of the junction and the bone reference point independent of tissue deformation and sensor tilting during contraction. Note the extent of dispersion of the myotendinous junction displacement with the current measurement setup (Figs 12A, 13A). The use of volume data could decrease the error in myotendinous junction displacement measurement.

Extraction of the myotendinous junction displacement parallel to the humerus

From the ultrasonic measurements, three features where extracted manually and with a (semi-) automatic approach. These features were the tracked movements of, first, the myotendinous junction, second, the humerus bone surface (bone line), and third, a reference point on the bone surface. The combination of these features allowed a calculation of the myotendinous junction movement in a bone local coordinate system, unlike other works in which the myotendinous junction was tracked in a skin local coordinate system (e.g. [9, 10] used ultrasonic reflective markers on the skin as reference points).

Of these three features, the automatic tracking methods for myotendinous junction and bone orientation via bone line tracking both proved to be feasible. It should be noted, that during preprocessing of the ultrasonic images a Gaussian filter was used, which could potentially contribute to some amount of error in the detection of the bone line and bone reference point. However, the automatic extraction of the bone reference point (represented by “some kind of texture” of the bone surface) based on the Lucas-Kanade was not used due to the correlation of the difference between manual and Lucas-Kanade tracking with the wrist force (R = −0.58). As the investigator selected distinct reference points (e.g. where brachialis fibers connect with the bone), the manual labeling was assumed to be the ground truth. This assumption led to the conclusion that there was a systematic error proportional to the wrist force F_w in case of Lucas-Kanade tracking of the bone reference point.

As a consequence, the tracking of the bone reference point, as part of the procedure to estimate , was performed manually according to Eq (31). Due to this one manual step, the former automatic method was regarded as (semi-) automatic.

In this study the myotendionous junction is observed at the point where the aponeurosis forks into the muscle (here described as myotendionous junction). The movement of the myotendionous junction is than solely attributed to changes in tendon length during model building. Further work could include more complex model configurations of the myotendinous junction such as the aponeurosis and fascia.

As mentioned before, the usage of a handheld linear ultrasonic probe limited the performance of data extraction for the myotendinous junction. It seems plausible, that this limitation also affected the tracking of the bone reference point. The bone was the structure with the greatest distance to the ultrasonic sensor. Therefore, tilting of the sensor probe had a larger effect on bone imaging as compared to the myotendinous junction. Also here, the usage of a two-dimensional ultrasonic sensor array to get volumetric data could be beneficial in future. The fitting of volumetric models of anatomical features could lead to a reliable feature tracking for all three anatomic features (bone orientation, bone reference point and myotendinous junction).

Tendon model and its fitting

During model building, three major assumptions were made that strongly influenced the observed model behaviour. First, there is no preloading of the tendon. Second, the effective lever length l_eff is calculated based on a two-dimensional tendon path model as depicted in Fig 10. Third, the resulting elbow torque is solely generated by the biceps brachii (τ_r = 0).

The effect of l_eff miscalculation.

If the effective lever length of the biceps brachii was miscalculated this would lead to an elbow angle dependent over- or underestimation of the tendon force F_T1. A quantitative analysis was not applicable as l_eff could not be obtained from the participants in the presented measurement setup (see Fig 3). Direct in vitro measurements of the effective lever length of the lower biceps brachii tendon exist. Murray et al. presented data of the peak effective moment arm of the distal biceps brachii tendon over the length of the radius [42]. This data is depicted in Fig 14A. In addition, the peak effective moment arm of subject 0 to 2 over their radius length l_r was entered into the graph. The data of subject 0 to 2 lay within the trend of [42]. Therefore, it was assumed that the maximal moment arm (between 60^∘ and 100^∘, according to [42]) and effective moment arms for angles near the maximum were predicted reasonably well by the presented model.

The effect of τ_r on the measurement.

During model building, τ_r (as defined in Eq (36)) was assumed to be zero as the participant had been advised to flex the elbow joint with no co-contraction. But during the experiment it turned out that subjects had the tendency to co-contract to reach the desired force levels. If co-contraction was noticed during the experiment, it was pointed out to the subject with the request to stop co-contracting even if the requested force level could not be reached exactly. If other than the assumed muscles for co-contraction were dominant, the actual residual torque would be negative ().(46)

Due to Eq (46), co-contraction would lead to an underestimate of the torque generated by the biceps brachii () and therefore also an underestimate of Δl_mt. Analog an overestimate of the torque generated by the biceps brachii () is possible, if additional flexor muscles were engaged during contraction τ_r > 0. Whenever ΔL_mt was overestimated this seemed to be the dominant error source. During the fit process shrinking sub-ranges of Θ were fitted and led to increasing optimal fit results (in case of manual data extraction). Shrinking sub-ranges mean that lower angles were excluded from the fit. Considering that the biceps brachii has its approximate optimal working point at Θ ≈ 115^∘ [43], it seemed to be plausible that the biceps brachii was not the mainly engaged muscle for small elbow angles (arm extended).

The effect of tendon preloading on elbow torque estimation based on wrist force measurement.

The preloading of the tendon is not directly measurable in vivo. A deceptive way to estimate preloading would be to observe the force at the wrist during the unloaded state. However, even if the only physiological torque source was the biceps brachii the preloading would not be distinguishable from preloading of the wrist strap (cmp. Fig 3). If additional torque sources were assumed, the lower tendon of the biceps brachii could be preloaded even if the resulting torque sum at the elbow joint would be zero.

Therefor, tendon pre-loading remains to be unknown. If tendon pre-loading occurs, it would lead to a gradual shift towards the linear region (see Fig 2).

Discussion of the tendon force estimation.

However, during model building it was assumed that the biceps brachii was the only active muscle. This assumption might have led to an overestimation of the tendon force and, as a consequence, to an overestimation of the tendon stiffness. Shukla et al. [4] investigated the properties of the biceps brachii tendon in vitro. They found a mean failure load force of the tendon-radius insertion at 90^∘ elbow flexion of 360N (SE 120N). In this study, a mean tendon force during maximal voluntary contraction at 90^∘ elbow flexion of (SE 84N) was found. This force seems to be too high considering the failure load found by Shukla et al. [4]. However, the comparison with Shukla et al. is difficult because the tearing of the tendon from the bone was measured in vitro. In vivo, the distal biceps brachii tendon not only inserts into the bone but also into the forearm fasciae, which leads to a presumably higher failure force. In addition, other flexor muscles, such as brachialis and radio brachialis, may also contribute to an overestimation of the tendon force. However, the actual force generated by the biceps brachii could be estimated by using the ratio of the maximal cross sectional area of the biceps brachii and the combined maximal area of all elbow flexor muscles. This is based on the assumption that the cross sectional area of a muscle is directly proportional to the force generated during maximal voluntary contraction. According to Erskine et al. [44], the biceps brachii makes up 41.8% (SE 1.2%) of the overall maximal cross sectional areas of all flexor muscles. Following this line of thoughts, the biceps brachiii would exert a force of 580N (SE 50N). If it is assumed that this force was distributed to the bone and fascial insertions, a reasonable maximal voluntary force is reached which is comparable to the findings of Shukla et al.. Without such an estimation of the force generated by the biceps brachii, the tendon force is overestimated. An overestimate of the tendon force leads to an overestimate of the tendon stiffness.

In the captured ultrasonic videos it was noticeable that the brachialis muscle group changes its diameter. This, however, speaks in favour of an increasingly negative τ_r for smaller elbow angles. Therefore, the parameterized tendon model was not used to determine tendon stiffness (change in tendon force per change in tendon length). Rather, the presented model allowed the tendon length change to be determined based on the elbow torque. Further work should focus on the contribution of the remaining flexor muscles on the resulting elbow torque to predict τ_r more accurately.

Smart et al. have conducted a comparable study. In that study, the measured elbow torque was also attributed only to the contraction of the biceps brachii [9]. Their results should therefore be comparable to those in this work. They defined the tendon stiffness as the slope of a linear displacement over force function. This function was defined to include the tendon displacement at 20% and 80% of maximum voluntary contraction. For comparison, the same linear stiffness was calculated for the fit results of this study as demonstrated in Fig 14B and listed in 1. Smart et al. found a stiffness of for 9 young participants [9]. In the present study, for the optimal fit of each subject’s elbow angle subset (Table 1) a mean stiffness of () (in case of the automatic data extraction) and () (for the manually extracted data) was found. These linear stiffness results are considered not to be significantly different (two-sided T-test) as compared to the findings of Smart et al. [9] ( for the automatically extracted data, and for the manually extracted data). A general drawback of the linear stiffness even as medical parameter is, that small variances in 20% of can lead to large variances in c_lin,20-80 due to the non-linearity in the toe region (cmp. Fig 14B).

Conclusion

In the design of biomechanical models of a myotendinous apparatus for the prediction of joint angles during limb movements, an exact separation between muscle and tendon length is necessary. If the tendon length (i.e. tendon displacement from resting length) was estimated by means of a tendon model, the muscle length could be calculated. It is important that the tendon model includes a non linear toe region for small forces. Existing work used biomechanical data of tendon stiffness from anatomical cadaver or in vitro studies [1, 19, 21], whose transfer ability to living systems can be discussed in general. The present work proposed an experimental setup in which the biceps brachii distal myotendinous junction displacement is observed for different elbow angles at varying wrist forces (muscle contraction) in vivo. The myotendinous junction displacement in a bone local coordinate system was extracted from ultrasonic B-mode videos using the proposed manual and (semi-) automatic method. This myotendinous junction displacement was then used together with the wrist force and a two dimensional tendon path model to fit a non linear stiffness model motivated by the experiments of [2]. The two dimensional tendon path model showed good agreement with aspects of in vitro measurements (see Fig 14A). Due to the assumption, that the biceps brachii is the only active muscle, discussed in the previous section, the existing tendon model can not be used to extract the stiffness of the tendon. It is however suitable, for the prediction of biceps brachii tendon displacement at a given elbow torque for a given elbow angle sub-range. The fitted tendon model conformed with linear stiffness measurements by [10], who applied the same assumptions as in this work. Compared to [10], the presented tendon model is more suitable for the prediction of myotendinous junction displacement at small elbow torques than those models that only use a linear stiffness curve defined by two points (20% and 80% of the maximum voluntary force, [10]).

Supporting information

S1 Table. Anatomical measurements distances between anatomical landmarks.

https://doi.org/10.1371/journal.pone.0275128.s001

S2 Table. Bone line detection Canny-Edge-Detection parameters as in the OpenCV python package opencv-python version 4.5.2.54 at PyPI.

https://doi.org/10.1371/journal.pone.0275128.s002

S3 Table. Bone line detection Hough-Line-Transformation parameters as in the OpenCV python package opencv-python version 4.5.2.54 at PyPI.

https://doi.org/10.1371/journal.pone.0275128.s003

S4 Table. Bone line detection DBSCAN parameters as in the SciPy python package scipy version 1.6.3 at CONDA-FORGE.

https://doi.org/10.1371/journal.pone.0275128.s004

S5 Table. Fitted model parameters model parameters for the fits in Figs 12, 13 for each subject.

https://doi.org/10.1371/journal.pone.0275128.s005

S6 Table. Elbow angle dependent lever arms the distance between the elbow joint and wrist strap depends on the elbow angle (cmp. Eq (33)). Here this distance is listed for each subject.

https://doi.org/10.1371/journal.pone.0275128.s006

(PDF)

Citation: Mechtenberg M, Grimmelsmann N, Meyer HG, Schneider A (2022) Manual and semi-automatic determination of elbow angle-independent parameters for a model of the biceps brachii distal tendon based on ultrasonic imaging. PLoS ONE 17(10): e0275128. https://doi.org/10.1371/journal.pone.0275128

About the Authors:

Malte Mechtenberg

Roles: Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

E-mail: [email protected]

Affiliation: Biomechatronics and Embedded Systems Group, Bielefeld University of Applied Sciences, Bielefeld, NRW, Germany

https://orcid.org/0000-0002-8958-0931

Nils Grimmelsmann

Roles: Investigation, Methodology, Writing – review & editing

Affiliation: Biomechatronics and Embedded Systems Group, Bielefeld University of Applied Sciences, Bielefeld, NRW, Germany

Hanno Gerd Meyer

Roles: Data curation, Software, Writing – original draft, Writing – review & editing

Affiliation: Biomechatronics and Embedded Systems Group, Bielefeld University of Applied Sciences, Bielefeld, NRW, Germany

Axel Schneider

Roles: Conceptualization, Funding acquisition, Methodology, Project administration, Supervision, Writing – review & editing

Affiliation: Biomechatronics and Embedded Systems Group, Bielefeld University of Applied Sciences, Bielefeld, NRW, Germany

https://orcid.org/0000-0002-6632-3473

References

1. Zajac FE. Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Critical reviews in biomedical engineering. 1989;17(4):359–411. pmid:2676342

2. Rack P, Westbury D. Elastic properties of the cat soleus tendon and their functional importance. The Journal of Physiology. 1984;347(1):479–495. pmid:6231373

3. Duenwald S, Kobayashi H, Frisch K, Lakes R, Vanderby RJ. Ultrasound echo is related to stress and strain in tendon. Journal of biomechanics. 2011;44(3):424–429. pmid:21030024

4. Shukla DR, Morrey BF, Thoreson AR, An KN, O’Driscoll SW. Distal biceps tendon rupture: an in vitro study. Clinical biomechanics (Bristol, Avon). 2012;27(3):263–267. pmid:22030096

5. Idler CS, Montgomery WHr, Lindsey DP, Badua PA, Wynne GF, Yerby SA. Distal biceps tendon repair: a biomechanical comparison of intact tendon and 2 repair techniques. The American journal of sports medicine. 2006;34(6):968–974. pmid:16476918

6. Zhang Q, Adam NC, Hosseini Nasab SH, Taylor WR, Smith CR. Techniques for In Vivo Measurement of Ligament and Tendon Strain: A Review. Annals of Biomedical Engineering. 2021;49(1):7–28. pmid:33025317

7. Chaudhry S, Morrissey D, Woledge R, Bader D, Screen H. Eccentric and Concentric Exercise of the Triceps Surae: An In Vivo Study of Dynamic Muscle and Tendon Biomechanical Parameters. Journal of applied biomechanics. 2014;31(2):69–78. pmid:25322475

8. Lichtwark G. In: Müller B, Wolf SI, Brueggemann GP, Deng Z, McIntosh A, Miller F, et al., editors. Ultrasound Technology for Examining the Mechanics of the Muscle, Tendon, and Ligament. Springer International Publishing; 2017. p. 1–20.

9. Smart R, Baudry S, Fedorov A, Kuzyk S, Jakobi J. Influence of biceps brachii tendon mechanical properties on elbow flexor force steadiness in young and old males. Scandinavian Journal of Medicine & Science in Sports. 2017;28(3):983–991.

10. Smart R, Richardson C, Jakobi J. Repeatability of ultrasound in assessment of distal biceps brachii tendon. International Journal of Diagnostic Imaging. 2017;4(2).

11. Wiesinger HP, Kösters A, Müller E, Seynnes O. Effects of Increased Loading on In Vivo Tendon Properties: A Systematic Review. Medicine and science in sports and exercise. 2015;47(9):1885–1895. pmid:25563908

12. Kleiber T, Kunz L, Disselhorst-Klug C. Muscular coordination of biceps brachii and brachioradialis in elbow flexion with respect to hand position. Frontiers in Physiology. 2015;6. pmid:26300781

13. An KN, Hui F, Morrey BF, Linscheid RL, Chao EY. Muscles across the elbow joint: a biomechanical analysis. Journal of biomechanics. 1981;14(10):659–669. pmid:7334026

14. Cucca Y, McLay S, Okamoto T, Ecker J, Mcmenamin P. The biceps brachii muscle and its distal insertion: Observations of surgical and evolutionary relevance. Surgical and radiologic anatomy. 2010;32(4):371–375. pmid:19847376

15. Thorpe CT, Birch HL, Clegg PD, Screen HRC. Tendon Physiology and MechanicalBehavior: Structure–Function Relationships. In: TENDON REGENERATION Understanding TissuePhysiology and Developmentto Engineer Functional Substitutes. Elsevier Inc.; 2015. p. 3–39.

16. Bennett M, Ker R, Imery NJ, Alexander RM. Mechanical properties of various mammalian tendons. Journal of Zoology. 1986;209(4):537–548.

17. Hooley C, Cohen R. A model for the creep behaviour of tendon. International Journal of Biological Macromolecules. 1979;1(3):123–132.

18. Moriya T, Zhao C, An KN, Amadio PC. The effect of epitendinous suture technique on gliding resistance during cyclic motion after flexor tendon repair: a cadaveric study. The Journal of hand surgery. 2010;35(4):552–558. pmid:20189323

19. Buchanan TS, Lloyd DG, Manal K, Besier TF. Neuromusculoskeletal modeling: estimation of muscle forces and joint moments and movements from measurements of neural command. Journal of applied biomechanics. 2004;20(4):367–395. pmid:16467928

20. Proske U, Morgan D. Tendon stiffness: methods of measurement and significance for the control of movement. A review. Journal of biomechanics. 1987;20(1):75–82. pmid:3558432

21. Thelen DG. Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults. J Biomech Eng. 2003;125(1):70–77. pmid:12661198

22. Mechtenberg M, Grimmelsmann N, Meyer HG, Schneider A. 2-D ultrasound videos of the distal biceps brachii myotendinous junction displacement over varying elbow angles and tendon loads; 2022. Available from: https://doi.org/10.5281/zenodo.7081639.

23. Knudsen AB, Larsen M, Mackey AL, Hjort M, Hansen KK, Qvortrup K, et al. The human myotendinous junction: an ultrastructural and 3D analysis study. Scandinavian journal of medicine & science in sports. 2015;25(1):e116–e123. pmid:24716465

24. Leitner C, Jarolim R, Konrad A, Kruse A, Tilp M, Schröttner J, et al. Automatic Tracking of the Muscle Tendon Junction in Healthy and Impaired Subjects using Deep Learning. In: 2020 42nd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). IEEE; 2020. p. 4770–4774.

25. Krupenevich RL, Funk CJ, Franz JR. Automated analysis of medial gastrocnemius muscle-tendon junction displacements in heathy young adults during isolated contractions and walking using deep neural networks. Computer methods and programs in biomedicine. 2021;206:106120. pmid:33991901

26. Cenni F, Bar-On L, Monari D, Schless SH, Kalkman BM, Aertbeliën E, et al. Semi-automatic methods for tracking the medial gastrocnemius muscle–tendon junction using ultrasound: a validation study. Experimental physiology. 2020;105(1):120–131. pmid:31677311

27. Zhou GQ, Zhang Y, Wang RL, Zhou P, Zheng YP, Tarassova O, et al. Automatic myotendinous junction tracking in ultrasound images with phase-based segmentation. BioMed research international. 2018;2018. pmid:29750152

28. Lee SS, Lewis GS, Piazza SJ. An algorithm for automated analysis of ultrasound images to measure tendon excursion in vivo. Journal of applied biomechanics. 2008;24(1):75–82. pmid:18309186

29. Cronin NJ, Carty CP, Barrett RS, Lichtwark G. Automatic tracking of medial gastrocnemius fascicle length during human locomotion. Journal of applied physiology. 2011;111(5):1491–1496. pmid:21836045

30. Mathis A, Mamidanna P, Cury KM, Abe T, Murthy VN, Mathis MW, et al. DeepLabCut: markerless pose estimation of user-defined body parts with deep learning. Nature neuroscience. 2018;21(9):1281–1289. pmid:30127430

31. Sandler M, Howard A, Zhu M, Zhmoginov A, Chen LC. Mobilenetv2: Inverted residuals and linear bottlenecks. In: Proceedings of the IEEE conference on computer vision and pattern recognition; 2018. p. 4510–4520.

32. Deng J, Dong W, Socher R, Li LJ, Li K, Fei-Fei L. Imagenet: A large-scale hierarchical image database. In: 2009 IEEE conference on computer vision and pattern recognition. Ieee; 2009. p. 248–255.

33. Tomar S. Converting video formats with FFmpeg. Linux Journal. 2006;2006(146).

34. Canny J. A Computational Approach to Edge Detection. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1986;PAMI-8(6):679–698. pmid:21869365

35. Bradski G. The OpenCV Library. Dr Dobb’s Journal of Software Tools. 2000;.

36. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods. 2020;17(3):261–272. pmid:32015543

37. Ester M, Kriegel HP, Sander J, Xiaowei X. A density-based algorithm for discovering clusters in large spatial databases with noise. Proceedings of the Second International Conference on Knowledge Discovery and Data Mining. 1996;96:226–231.

38. Lucas BD, Kanade T. An Iterative Image Registration Technique with an Application to Stereo Vision. In: Proceedings of the 7th International Joint Conference on Artificial Intelligence—Volume 2. IJCAI’81. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc.; 1981. p. 674–679.

39. Shi J, Tomasi. Good features to track. In: 1994 Proceedings of IEEE Conference on Computer Vision and Pattern Recognition; 1994. p. 593–600.

40. Bouguet JY, et al. Pyramidal implementation of the affine lucas kanade feature tracker description of the algorithm. Intel corporation. 2001;5(1-10):4.

41. Sinelnikov RD. Atlas of Human Anatomy. mir publischer moscow; 1988.

42. Murray WM, Buchanan TS, Delp SL. Scaling of peak moment arms of elbow muscles with upper extremity bone dimensions. Journal of biomechanics. 2002;35(1):19–26. pmid:11747879

43. Chang YW, Su FC, Wu HW, An KN. Optimum length of muscle contraction. Clinical biomechanics (Bristol, Avon). 1999;14:537–542. pmid:10521638

44. Erskine RM, Fletcher G, Folland JP. The contribution of muscle hypertrophy to strength changes following resistance training. European journal of applied physiology. 2014;114(6):1239–1249. pmid:24610245