A lumped stiffness model of intermuscular and extramuscular myofascial pathways of force transmission

Mechanical behavior of skeletal muscles is commonly modeled under the assumption of mechanical independence between individual muscles within a muscle group. Epimuscular myofascial force transmission via the connective tissue network surrounding a muscle challenges this assumption as it alters the force distributed to the tendons of individual muscles. This study aimed to derive a lumped estimate of stiffness of the intermuscular and extramuscular connective tissues and to assess changes in such stiffness in response to a manipulation of the interface between adjacent muscles. Based on in situ measurements of force transmission in the rat plantar flexors, before and after resection of their connective tissue network, a nonlinear estimate of epimuscular myofascial stiffness was quantified and included in a multi-muscle model with lumped parameters which allows for force transmission depending on the relative position between the muscles in the group. Such stiffness estimate was assessed for a group with normal intermuscular connective tissues and for a group with increased connectivity, mimicking scar tissue development. The model was able to successfully predict the amount of epimuscular force transmission for different experimental conditions than those used to obtain the model parameters. The proposed nonlinear stiffness estimates of epimuscular pathways could be integrated in larger musculoskeletal models, to provide more accurate predictions of force when effects of mechanical interaction or altered epimuscular connections, e.g. after surgery or injury, are substantial.


Introduction
It is traditionally thought that skeletal muscles transmit force to their insertions on bones primarily through the myotendinous junction (Tidball 1991). This common understanding has been challenged by the evidence of secondary pathways of force transmission (Huijing 2009;Maas and Sandercock 2010). According to this new perspective, part of the force produced by a muscle can be exerted on neighboring muscles and surrounding extramuscular structures through a continuous network of connective tissues. Such a network extends from the endomysial-perimysial layers (Trotte et al. 1995;Purslow 2010) to the epimysium of muscles and further to the connective tissues between adjacent muscle bellies (intermuscular) and tissues supporting nerves and blood vessels (extramuscular) (Maas and Sandercock 2010;Higham and Biewener 2011). These inter-and extramuscular pathways provide the means for epimuscular myofascial force transmission (Huijing et al. 1998;Maas et al. 2001;Yucesoy et al. 2006), in addition to the parallel myotendinous path to bone. As a consequence, muscles do not necessarily function as independent actuators.
Until recently, multi-muscle models considered skeletal muscles as independent actuators (Woittiez et al. 1985;Zajac and Winters 1990;Jacobs et al. 1996;Raasch et al. 1997;Correa et al. 2011;Lee et al. 2015). Moreover, these models have commonly been parameterized by dissecting muscles free from the structures surrounding them (Close 1964;Rack and Westbury 1969;Burke et al. 1976;Woittiez et al. 1985;Ettema et al. 1990). Therefore, stiffness of inter-or extramuscular pathways is neglected. A first formal description of extramuscular linkages was encountered in a phenomenological Hill-type muscle model to explain an unexpected, and first ever reported, difference between forces measured at both ends of a muscle (Devasahayam and Sandercock 1992). Despite its novelty, the authors considered this model speculative, as it was not quantitatively tested nor was a stiffness estimate for extramuscular linkages reported. More recently, a specific Hill-type model was extended to approximate the gearing between transverse and longitudinal forces on the muscle belly. Although this model was able to predict changes in muscle force produced by external loads exerted on the muscle belly, the muscle preparation and the resulting parametric model of transverse loading used a single muscle isolated from its surrounding (Siebert et al. 2014a, b). In a preliminary attempt to illustrate, although only qualitatively, the effects of extraand intermuscular connections, a physical model representing connections between rat ankle dorsi-flexors was devised using springs elements of different stiffness and inextensible linking elements (Maas et al. 2001). Finally, finite element (FE) modeling emerged as a promising approach to account for several passive properties of muscle tissue as a continuum. These models allowed to quantitatively assess mechanisms of force transmission and the effects of mechanical interaction between muscle fibers at the extracellular matrix interface on sarcomere length, shear strain and force (Yucesoy et al. 2002;Sharafi and Blemker 2011;Zhang and Gao 2012), and between whole muscles linked by intermuscular elastic structures (Yucesoy et al. 2003). However, the computational load of FE models makes the integration of such multi-muscle systems into whole limb or body musculoskeletal models more difficult. A simpler model accounting for macroscopic muscle behavior seems, at first instance, preferable. Moreover, the complexity of intermuscular connections and macroscopic differences between muscle groups, e.g. moment arm, number of joints spanned, distance between origin-insertion of neighboring muscles, presents an open challenge for computing local stress-strain maps of the synergistic muscle group as a whole. Also, previous studies have shown that the extent of epimuscular myofascial effects can vary between muscle groups (see Table 1 in Huijing 2009) and species (Maas and Sandercock 2008), which requires modeling tools to be adaptable. Despite the fact that a number of studies have described force transmission via epimuscular myofascial linkages, no stiffness values have been assessed directly from experimental data. Note that the stiffness of epimuscular linkages in above described FE models was the result of optimization (Yucesoy et al. 2002;Maas et al. 2003).
We have shown that the ratio of transmitted force between myotendinous and epimuscular (inter-or extramuscular) pathways changes due to altered mechanical properties of connective tissues between muscles (Bernabei et al. 2016). Such changes can be a consequence of muscle injury (Garrett et al. 1988;Kääriäinen et al. 2000;Silder et al. 2010), disease (Booth 2001;Smith et al. 2011), surgical treatment (Almeida-Silveira et al. 2000Smeulders and Kreulen 2007) or aging (Kovanen et al. 1987;Ramaswamy et al. 2011). To take into account effects of post-operative scar tissue development after tendon transfer of spastic muscles, a secondary non-myotendinous pathway of force transmission was added to an EMG-driven musculoskeletal model of the quadriceps muscles (Andersen 2009). Although this model provides a simple solution to integrate force transmission in pathological conditions, the ratio of force transmitted to the new insertion and via scar tissue to the patella was arbitrarily set, i.e., not determined by the mechanical properties of the post-operative scar tissues.
The aims of this study were (i) to derive a lumped estimate of stiffness of the intermuscular and extramuscular connective tissues and (ii) to assess changes in such stiffness in response to a manipulation of the interface between adjacent muscles. Based on in situ measurements of force transmission in the rat soleus (SO) and lateral gastrocnemius and plantaris complex (LG + PL), before and after resection of their connective tissues network, a nonlinear estimate of stiffness was quantified and included in a phenomenological multi-muscle model with lumped parameters. We used this model to predict the ratio of force transmitted via epimuscular myofascial and myotendinous pathways of SO and LG + PL for both normal and enhanced intermuscular connectivity. Finally, the relative contribution of extramuscular pathways of force transmission, i.e., SO and LG + PL neurovascular tracts, with respect to the epimuscular myofascial and myotendinous pathways was assessed.

Animals
Experiments were performed on male Wistar rats (n = 14, mean ± SD body mass 308.0 ± 11.1 g), which were divided into two groups: the first group had normal connective tissues (n = 7, normal connectivity group, NO), while the second one was tested after manipulation of the connectivity between SO and LG + PL (n = 7, tissue-integrating mesh group, TI). All surgical and experimental procedures were approved by the Committee on the Ethics of Animal Experimentation at the VU University Amsterdam and in strict agreement with the guidelines and regulations concerning animal welfare and experimentation set forth by Dutch law.

Surgical procedures for intermuscular connectivity manipulation
One to two weeks prior to the in situ measurements, TI group animals were implanted with a tissue-integrating surgical mesh at the muscle belly interface of SO and LG + PL. Surgical procedures for the implantation have been previously reported in detail (Bernabei et al. 2016). Briefly, following preparation for aseptic surgery and inhalation anesthesia (2-3 % isoflurane), the lateral side of the posterior crural compartment was accessed by a partial fasciotomy of the surrounding crural fascia, limiting the exposure of the posterior crural muscles to the proximal two-thirds of the SO and LG muscle-tendon units (MTUs). The distal portion of the biceps (accessory head) and posterior crural fascia, covering the distal myotendinous junctions of SO, LG, medial gastrocnemius and PL as well as the Achilles tendon, were left intact. After blunt dissecting the connective tissues between the dorsal side of the SO and the ventral side of LG + PL, a tissue-integrating mesh (n = 7; Premilene mesh, B. Braun Melsungen AG, Germany) was sutured to the SO muscle belly at the interface of SO and LG + PL muscles. Meticulous dissection prevented damage of the SO neurovascular tract (NV1), which runs centrally between SO and LG muscle bellies. A one-time pre-operative subcutaneous injection of a pain-killer (0.02 mg/kg, Temgesic , Schering-Plough, Maarssen, The Netherlands) was administered to prevent discomfort of the animal following the survival surgery. Additional doses were given 1-2 days after the surgery if signs of pain were noticed.

Surgical procedures for in situ measurements
According to standard procedures in our laboratory (Maas et al. 2001), the animals were anesthetized by an intraperitoneal injection of urethane solution (1.2 ml/100 g body mass, 12.5 % urethane solution). Extra doses were administered if necessary (1.5 ml maximally) until withdrawal reflexes to a pain stimulus were suppressed completely. At the end of the measurements, the animals were euthanized with an overdose of pentobarbital sodium (Euthasol 20 %) injected intracardially, followed by double-sided pneumothorax. Surgical procedures and technical specifications of the experimental setup for the in situ measurements have been described in detail elsewhere (Bernabei et al. 2015). The posterior crural compartment was exposed and biceps femoris and medial gastrocnemius muscles were removed. SO, LG and PL muscles were, as a group, dissected free from surrounding structures, preserving the bone insertions as well as the myofascial connections at the interface between the SO and LG + PL muscle bellies (Fig. 1a). The neurovascular tracts supplying this muscle group were left intact. The distal tendon of SO was dissected free from the rest of the Achilles tendon. Kevlar wires were used to connect the proximal and distal tendons of LG + PL as well as the distal tendon of SO to force transducers, which were positioned in such a way that forces could be measured in the muscle's line of pull.
It should be noted that the SO origin was left intact to preserve the original muscle position relative to the tibia; thus, SO proximal tendon force was not measured. A bipolar cuff electrode connected to a constant current source was folded around the sciatic nerve. The common peroneal nerve was severed, leaving only the tibial nerve and the sural branch intact.

Nerve stimulation
SO and LG + PL were excited simultaneously via supramaximal stimulation of the sciatic nerve (0.4 ± 0.1 mA, 500 ms, 100 Hz). Two twitches were evoked 1.5 and 1 s before each tetanic stimulation. Muscles were allowed to recover for 2 min between subsequent stimulations.

Experimental protocol
The rat was mounted in the experimental setup by clamping the femur and the foot such that knee and ankle joints were kept at 90 • , which was defined as the reference position (P REF ). Markers were placed near the SO distal and LG + PL proximal and distal tendons to apply different relative positions between SO and LG + PL MTUs. To prevent history effects, multiple contractions at high and low lengths were elicited prior to data collection. Isometric forces exerted at all three tendons were measured simultaneously for different positions of LG + PL relative to SO muscle (P R ). The proximal and distal LG + PL tendons were repositioned in steps of 1 mm from P REF − 3 mm in distal direction to P REF +3 mm in proximal direction (Fig 1b).
LG+PL tendons were repositioned together, while the distal and proximal tendons of SO were kept at the same position corresponding to P REF . Thus, MTU lengths were constant over all the imposed conditions, with SO and LG + PL operating on the ascending region of their length-force curves (mean ± SD at P REF : Bernabei et al. 2015). Following the same experimental protocol, tendon forces were measured for three different conditions to estimate the individual contribution of inter and extramuscular pathways: (i) intact conditions (NO, n = 4; TI, n = 7); (ii) resection of proximal and distal intermuscular myofascial connections (NO, n = 4; TI, n = 7); (iii) complete isolation of synergists by resecting the neurovascular tract between SO and LG + PL (n = 1). It should be noted that the latter resection was possible only on rats with unaltered connectivity, as in the TI group the blood supply and nerve to SO could be no longer distinguished from the bundle of newly grown intermuscular connective tissue. Moreover, this procedure caused the SO to be no longer excitable. In a subset of animals, only distal myofascial connections were resected (NO, n = 3) to assess the changes in tendon force when only the proximal, instead of proximal and distal, connective tissues pathway was involved (Fig. 1a).

Analysis of experimental data
Isometric forces were assessed from the force-time series: passive force was assessed by calculating the mean over a 50 ms time window before the tetanic stimulation and total force was assessed by calculating the mean over the last 50 ms of the tetanic stimulation (Fig. 2). Two estimates of intermuscular mechanical interaction were assessed: with repositioning both proximal and distal LG + PL tendons simultaneously, we calculated (i) changes in force exerted at the distal tendon of SO ( F SO ), with SO kept at a constant length and (ii) changes in the force difference between the proximal and distal tendons of LG + PL, with LG + PL kept at a constant length. Changes in SO force were expressed relative to the value at P REF ( The LG + PL force difference ( F LGPL ) was calculated by subtracting the force exerted distally from the force exerted LG + PL. Passive and total forces at soleus distal (SO), LG + PL proximal (LG + PL prox) and LG + PL distal (LG + PL dist) tendons were measured simultaneously before (passive) and during (active) tetanic contraction. Gray areas indicate the 50 ms time-windows for extraction of forces during passive and active state of muscles. Mechanical interaction between SO and LG + PL is noticeable in the first portion of SO relaxation, which is much quicker than what can be expected from the slow SO muscle fibers and, thus, attributed to relaxation of the faster LG + PL muscle fibers proximally (i.e., F PROXIMAL − F DISTAL ), so that a positive difference indicates a higher proximal force. This force difference is a direct measure of the magnitude of net epimuscular myofascial force transmission Maas and Sandercock 2010). F LGPL and F SO as defined above will be referred to as the calibration dataset.
To test for effects of muscle relative position on proximal and distal LG + PL forces, two-way analysis of variance (ANOVA) for repeated measures (within factor: muscle relative displacement) was performed. A linear regression analysis was used to evaluate the relation between F SO and F LGPL over P R ∈ [−3; +3] mm. F SO and F LGPL were fitted with a polynomial function, with the order determined by stepwise regression. Generalized estimating equations (GEE) were used to test for significant differences between groups including the interactions with the factor groups (NO, TI), to test for differences in the slope of the models. p values <.05 were considered statistically significant.

Inter-and extramuscular pathways of force transmission
The myofascial connective tissues at the interface between SO and LG + PL muscle bellies, classified as the intermuscular pathway of force transmission (INT), were preserved and represented in the model as nonlinear elastic elements (Fig. 1c). This pathway was divided into a proximal and a distal part with respect to the neurovascular tract, which runs centrally between SO and LG muscle bellies (see Bernabei et al. 2016 for a detailed description). Bundles of connective tissues surrounding blood vessels and nerves reaching the posterior crural muscles can be divided in two main pathways for extramuscular force transmission (Fig. 1a). The first pathway (NV1) runs from the proximal-ventral side of LG to the dorsal side of SO. The second pathway (NV2) contains branches of the popliteal artery and tibial nerve running from the popliteal fossa to the proximal-dorsal surface of lateral and medial gastrocnemius muscle bellies. A clear depiction of the tibial nerve branches innervating the superficial muscles of the posterior crural compartment has been reported previously (Tijs et al. 2014). Although NV1 was defined here as an extramuscular pathway of force transmission, because it ultimately leads forces outside the synergistic muscle group, it should be noted that it also contributes to intermuscular force transmission, since it establishes a direct connection between SO and LG epimysial layers.

Model description
A phenomenological lumped-parameter model was used to represent SO and LG + PL interactions (Fig. 1c). Each mus-cle was described by a contractile element (CE) surrounded, both in series (SE) and in parallel (PE), by connective tissue structures. In particular, the mechanical properties of (i) intermuscular connective tissues INT, i.e., proximal and distal epimuscular myofascial connections and (ii) extramuscular connective tissue bundles NV1 and NV2, were summarized by elastic elements with piecewise linear stiffness properties. A formal explanation of the equations and parameters used to describe the mechanical lumped-parameter system can be found in the "Appendix". Briefly, the non-zero force balance at the tendons of a non-isolated muscle actuator can be seen as a direct measure of the amount of force transmitted interand extramuscularly to surrounding structures (Maas et al. 2004). Thus, F LGPL was considered as an estimate of the force redirected from LG + PL to neighboring structures via all the considered inter-and extramuscular connections. It should be noted that no direct information was available on the actual force vectors associated with the inter-and extramuscular pathways. F proximal and F distal were measured and included in the model for the three different experimental conditions applied: (i) intact connective tissues pathways (intact), (ii) after resection of INT (post I) and (iii) after resection of INT and NV1 (post II). With intact intermuscular and extramuscular connective tissue pathways, the balance of force produced by active muscles within the system shown in Fig. 1c is given by: where F is the force redirected to the given node via the superscripted element. Because MTUs of LG + PL and SO were kept at the same length between conditions and muscle contractions were isometric, F CE and F PE can be assumed to be constant and the forces in the system can be resolved with a static force equilibrium. Following progressive resection of connective tissues pathways, this system of equations was modified by setting the amount of force running through resected pathways to zero. It should be noted that proximal and distal INT pathways are represented in Eq. 2 by a single variable F INT . Solving Eqs. 1-2 for F INT allowed to estimate the amount of force transmitted via intermuscular connective tissues pathways. A similar approach has been used previously to quantify myofascial force transmission between the multiple heads of the extensor digitorum longus muscle (Zhang and Gao 2014) and changes in the mechanical effect of flexor carpi ulnaris after dissection from surrounding connective tissues in the rat (Maas et al. 2013).

Stiffness estimates of intermuscular and extramuscular pathways
The stiffness estimate for each pathway of force transmission was based on input-output data obtained from the in situ experiments (calibration dataset). The mechanical properties of passive connective tissues, from tendon to blood vessels, can be characterized as quasi-static (Fung 1981), so that continuous force transmission changes occurring during dynamic tasks can be approximated by changes in isometric force measured at discrete relative positions. The behavior of collagen fibers, of different initial orientation, starting to stretch at different extensions can be represented with a piecewise linear function (Winters and Woo 1990). Therefore in this study, the estimate of stiffness (K ) for each interand extramuscular pathway of force transmission (INT, NV1, NV2) was defined as the ratio of force increase ( F) over a unit of relative displacement ( P R ). Based on the system of equations which describes force equilibrium for different conditions of connective tissues resection, estimates of the force transmitted for each SO/LG + PL relative displacement were calculated as: LGPL * (c−1) where F LGPL is the force difference at LG + PL proximaldistal tendons for the intact, post I and post II conditions. The dimensionless parameter c was set to 0.9 to describe the ratio between the length of NV1 before and after resection of myofascial linkages (see "Appendix"). Sensitivity analysis with c ∈ (0.05; 1) resulted in a root-mean-square error corresponding to maximally 6.0 % of the average force transmitted via INT pathways over the imposed range of muscle displacements. To describe the distribution of force transmitted via non-myotendinous and myotendinous pathways, ratios of force transmitted via the intermuscular and extramuscular pathways and mean LG + PL proximal tendon force were calculated (

Model validation
The lumped-parameter model was validated by comparing the predicted LG + PL tendon force balance ( F mod LGPL ) with previously collected experimental data (Bernabei et al. 2016), here referred as the testing dataset ( F exp LGPL ). The experimental protocol used for measuring forces in the testing dataset differed from that used for the calibration dataset. For the testing dataset, LG + PL force changes were measured applying a physiological range of lengths and relative displacements (P prox R ) between SO and LG + PL, with only LG + PL proximal tendon repositioned from P prox R − 3 mm to P prox R + 3 mm, while the distal tendons were kept at the same position. For the calibration dataset instead, both proximal and distal tendons of LG + PL were repositioned together from P R −3 mm to P R +3 mm, with constant muscle lengths for LG + PL and SO (see above). Predicted values of LG + PL proximal-distal force difference were calculated as follows: where P REF corresponded to P R = 0; this SO/LG + PL relative position represented the matching condition for muscle lengths and relative displacement between the calibration dataset and the testing dataset. P prox R defined the proximal displacement of LG + PL. The scale-factor r described the ratio of myofascially transmitted force when a proximal displacement only (testing dataset), rather than a proximal and distal displacement (calibration dataset), was applied.

Analysis of model results
ANOVA for repeated measurements (within factor: LG + PL proximal tendon position, between factor: predicted vs. measured data) was used to test for differences between the values of LG+PL force difference predicted by the model and those of the testing dataset. A Bland-Altman plot was used to evaluate the agreement between predicted and measured values of LG + PL force differences. The bias was calculated as the mean difference between measured ( F exp LGPL ) and predicted force values ( F mod LGPL ).

Tetanic forces
A linear relationship between F LGPL and F SO was found (Fig. 3a). For both experimental groups, regression analysis revealed a positive slope (model fit: p < .001) with a high degree of correlation (NO: R = 0.98; TI: R = 0.95). The positive linear regression between SO and LG + PL total forces suggests that the force discrepancy at LG+PL tendons was balanced by the force exerted at SO distal tendon, which can be explained by force transmitted from LG + PL to SO via connective tissues at their muscle belly interface.
GEE showed that the slope of the regression line for total force was significantly steeper in the TI group than in the NO group ( p = 0.001), suggesting that the ratio between the force transmitted via epimuscular pathways with respect to the myotendinous pathway was increased after the TI manipulation surgery. Compared to the NO group, the increased force range with the same applied relative displacement indicates that the stiffness of intermuscular linkages was enhanced in the TI group (NO: F LGPL = 0.51N, TI: F LGPL = 2.12N at P R = +3 mm). Finally, it should be noted that F LGPL reached 0 N at P R = 0mm (REF position) and became negative for P R ∈ [−3; 0] mm, which suggests a reversal of the sign of the epimuscular force transmission vector at P R = 0 mm.

Passive forces
The relationship between passive F LGPL and normalized SO force was clearly not linear for both NO and TI groups (Fig. 3b). This result suggests that in passive muscle conditions, the force difference at LG + PL proximal and distal tendons was balanced only to some extent by force exerted at the SO distal tendon, especially for P R < 0 mm. Such non-exact force matching can be explained by force transmission via the proximal SO tendon and other non-muscular structures, which were not measured. As we considered the correspondence of F LGPL with the normalized SO force an a priori condition for modeling force transmission via intermuscular connective tissues, estimates of epimuscular stiffness and tendon force predictions were not extended to the passive muscle state.

Estimates of force transmission via intermuscular and extramuscular pathways
Force generated by LG + PL muscle fibers was partially redirected to SO (via INT and NV1) and to non-muscular extramuscular structures (via NV2) with a variable extent as a function of relative position (Fig. 4). The amount of epimuscular myofascial force relative to myotendinous force (F/F prox LGPL ) for NO (Fig. 4a) and TI (Fig. 4b) was maximal at P R = −3 mm (7.7 %) and at P R = +3 mm (17.7 %), respectively. NV1 and NV2 contributed to transmit force up to 1.1 and −4.0 % of proximal F LGPL at P R = −3 mm, respectively. It should be noted that the force transmission vector via NV1 and NV2 was oriented in opposite direction for most conditions. An overview of the stiffness estimates derived from the calibration dataset is reported in Table 1. These results show that the force distribution among SO and LG + PL synergists is the result of a complex interaction between pathways of force transmission. The behavior of these pathways differed in terms of stiffness and configura- Fig. 3 Comparison of SO and LG + PL force discrepancies in active and passive state of muscles. Regression analysis between F SO (mean + SD) and F LGPL (mean + SD) with SO/LG + PL relative displacement from P R = −3 to P R = +3 mm. Total forces (a) and passive forces (b) are compared between the normal group (NO) and the tissue-integrating mesh group (TI). Linear regression lines are illustrated for each group (NO: continuous line; TI: dashed line). R-square values of the linear regression for the total force of NO and TI groups were 0.98 and 0.95, respectively. Passive forces in NO and TI were clearly nonlinear tion, i.e., angle and direction of the force transmission vector in response to changes in SO/LG + PL relative position.

Validation study
ANOVA showed no significant differences between F exp LGPL and F mod LGPL ( p = .956). In the testing dataset, LG + PL proximal-distal force difference varied with LG + PL proximal position P prox R for both NO and TI groups ( p < .001, Fig. 5a, b, respectively); this force variation did not significantly differ between predicted and measured values ( p = .167). The Bland-Altman plot shows a bias of the predicted dataset with respect to the measured one (NO: 0.106N, TI: −0.071N, or 13.1 and 8.8 % of the confidence interval, respectively). In particular, force values predicted by the model were overestimated at a low values of F LGPL for NO and underestimated at high F LGPL for TI (Fig. 5c). The limits of agreement between model and experimental data were comprised within a range of about 0.840 N (±1.96* SD = 0.475; −0.366 N).

Discussion
We here presented a phenomenological multi-muscle model, which takes into account the intermuscular and extramuscular pathways of force transmission of the rat plantar-flexor muscles. Nonlinear estimates of stiffness for these pathways were derived based on in situ measurements of force transmission. We showed that such estimates can be used to predict the amount of force transmitted via non-myotendinous pathways at physiological lengths and relative positions of mus-cles. Although the magnitude of non-myotendinous force transmission effects in in vivo conditions are disputed (Herbert et al. 2008;Maas and Sandercock 2008;Tijs et al. 2015b), imaging studies in humans do imply relevant effects  (N V 1, N V 2), over each step of SO/LG + PL MTUs displacement ( P R ). Force values are expressed as a percentage relative to the LG + PL proximal tendon total force ( F/F LG prox [%]) for the normal group (a) and the tissue-integration group (b). As an example, the distribution of force among the different pathways in respect to the myotendinous one, i.e., the proximal LG + PL tendon, is shown for P R = 2 mm (Bojsen-Moller et al. 2010;Huijing et al. 2011;Tian et al. 2012;Yaman et al. 2013). Epimuscular connective tissues pathways may substantially affect muscle mechanics when these connections are either stiffer or more compliant than normal, e.g. because of surgery, injury or pathology (Yucesoy et al. 2007;Smeulders and Kreulen 2007;Huijing et al. 2010;Maas and Huijing 2012b).
To our knowledge, this is the first model that relates the stiffness of connective tissues to changes in tendon forces, depending on the relative position of synergistic muscles. A difference between forces exerted simultaneously at the origin and insertion of a muscle is considered a direct measure of the force transmitted by a muscle to surrounding structures via non-myotendinous pathways Maas et al. 2001). The relative position of a muscle with respect to surrounding muscles affects muscle function by altering the amount of force transmitted epimuscularly Maas et al. 2004). In this study, we combined these features to obtain estimates of stiffness of intermuscular myofascial linkages and of neurovascular tracts. By deriving such estimates in a condition of normal (NO) and enhanced connectivity (TI), the extent of transmitted force could be related to the mechanical properties of these connective tissues. Moreover, given the correspondence between the imposed SO and LG + PL relative positions with in vivo knee and ankle angles (Johnson et al. 2011), changes in force transmission can be represented as a function of joint angle (Fig. 6a). By linearly scaling the stiffness estimates between normal and enhanced connectivity conditions, the resulting three-dimensional function (Fig. 6b) allows a novel and concise representation of the epimuscular K N : parameters describing the stiffness of a given force pathway over a unit of relative displacement as defined in Eq. 7

Fig. 5
Model predictions versus measured experimental data. Validation of the epimuscular force transmission model (n = 12, calibration dataset) using data from a previously reported in situ experiment as a testing dataset (n = 24) in which physiological lengths and relative position of SO and LG + PL were applied (Bernabei et al. 2016).
Averaged predicted values (black circles) of LG + PL proximal-distal force difference with repositioning LG + PL proximal tendon ( P prox R ) from −3 to +3 mm from P REF are reported together with the corresponding experimental values (red squares, mean ± SD) for animals with enhanced connectivity (a, n = 7) and animals with normal con-nectivity (b, n = 8). Bland-Altman plot (c) showing the average bias (NO: 0.105N, TI: −0.071N) of the modeled data. Each point represents a force value estimated for all animals within NO (circles) and TI (squares) groups. Values of F LG+PL on the x-axis were expressed as a percentage of the total range of LG + PL proximal-distal force difference (max( F LG+PL ) − min( F LG+PL )) to reduce variation between animals. The limits of agreement, expressed as the mean bias ± 1.96 SD, for individual force data are shown (lower limit: −0.475; upper limit: 0.366N) force transmission phenomena, which could be integrated in a more comprehensive musculoskeletal model. This would enable predictions of force changes at the tendons of nonisolated muscles in functional tasks when connective tissues are altered.

Limitations
Given the novelty of considering non-myotendinous pathways of force transmission in musculoskeletal modeling, a number of assumptions had to be made to reduce the LGPL ) plotted a as a function of knee angle corresponding to the applied proximal LG + PL length changes (P prox R ) and b as a function of both knee angle and connective tissues stiffness. In the top panel, predicted values of intermuscular and extramuscular force transmission for a 90 • knee angle are highlighted for the evaluated intermuscular connectivity conditions, i.e., normal (NO) and enhanced connectivity (TI). The secondary axis expresses the predicted F mod LGPL as a percentage of the mean LG + PL proximal tendon force (F LG+PL , proximal: 11.5 ± 0.2N). In the bottom panel, linear interpolation of force transmission estimates between NO (green) and TI (red) over the range of corresponding knee angles complexity of the multi-muscle system and to isolate the effects of inter-and extramuscular connective tissues. First, unlike previous modeling studies (Maas et al. 2003;Yucesoy et al. 2007), we assumed the consequences of changes in the distribution of lengths of sarcomeres with changes of muscle relative position on muscle force to be negligible, an assumption that has been verified experimentally for the muscles considered here (Tijs et al. 2015a). Second, although local deformations of muscle bellies and differential strain of all force transmission pathways, from muscle fibers to tendons and epimuscular connective tissues, are expected (Finni et al. 2003;Blemker et al. 2005;Chi et al. 2010;Tijs et al. 2015b), we assumed that muscle relative position was the sole determinant of force changes in the muscle group. In fact, these local effects occur as a consequence of changes in the imposed relative positions, so that relative position was selected as a global descriptor of length changes for the whole mechanical chain. It is conceivable that changes of relative position occurring at both tendons simultaneously would elicit tension of non-myotendinous linkages to a higher extent and more homogeneously than considering proximal or distal positional changes only. In this regard, a comparison between the forces measured at LG + PL tendon after resection of all intermuscular connective tissues and after resection of the distal intermuscular connective tissues only was included in the "Appendix" (Fig. 8). Third, the contribution of the Achilles tendon to the distribution of force within the muscle group (Sandercock 2000;Tijs et al. 2014) could not be taken into account, as the distal tendons of SO and lateral gastrocnemius were separated. This was necessary to measure individual tendon forces and to distinguish effects of non-myotendinous force pathways from the myotendinous ones. Finally, stiffness estimates were not representative of a passive muscle state, since this condition did not guarantee negligible deformation of muscle bellies during contraction. This was a necessary assumption of this model (see "Appendix"). In fact, it has been reported that changes in connective tissue layers between synergistic muscles are not equally expressed in passive and active states (Smeulders et al. 2002;Maas and Huijing 2012a;Bernabei et al. 2016). In agreement with these studies, we did not find equilibrium between F LGPL and F SO in the passive condition (Fig. 2b). Given the aforementioned assumptions and limitations, tendon forces predicted by this model should be carefully interpreted, especially when considering a different anatomy for connective tissues connections and various muscle activation levels.

Implications for musculoskeletal modeling
When the assumption of mechanical independence of muscles cannot be verified, the output of a model based on this assumption may fail to predict forces in the musculoskeletal system represented. This has important implications for the experimental techniques used to model length-force characteristics of human muscles in vivo (Herzog and Keurs 1988;Hoang et al. 2007;Winter and Challis 2010). When using moment-angle data to reconstruct length-force curves of the muscles spanning a given joint, changes in relative position of a one-joint muscle respect to a two-joint muscle and the location of force measurement, i.e., proximal or distal tendons, can affect the result (Bernabei et al. 2015). The approach described in this study allows to explain inconsistent or variable outcomes of in vivo experiments when the assumption of mechanical independence is not guaranteed.
Considering in situ experiment instead, very often these are performed on non-isolated muscles or muscle compartments as well as dissected muscles still connected to their blood supply and innervation (Burke et al. 1976;Scott et al. 1996;Whitehead et al. 2001). These preparations are prone to effects of inter-and extramuscular force transmission, which may mask the actual mechanical behavior of the muscle fibers or introduce a force offset in between conditions with different relative positions between muscles. Given the results of the present study, a preliminary test on force differences at the tendons of a muscle with imposing different lengths or relative positions may help assessing the extent of myofascial effects. Alternatively, we showed that passive elastic elements are sufficient to model muscles' interaction when the assumption of mechanical independence cannot be verified. Such a conclusion is in agreement with the study of Devasahayam and Sandercock (1992), where extramuscular passive elements were introduced in a simple Hill-type model with isolated actuators to explain force discrepancies at both ends of a muscle.
FE multi-muscle models have been used to represent effects of extramuscular force transmission on muscle mechanics with changes of muscles relative position (Maas et al. 2003;Yucesoy et al. 2003). However, the FE method is unnecessarily complex for describing global tendon forces, and too computationally demanding to be embedded in largescale biomechanical models of the musculoskeletal system. In contrast, the model described here can be used to predict force changes at the muscle and muscle group level, thus matching the output of musculoskeletal models representing muscle function at a macroscopic level.

Implications for modeling of in vivo human function and pathology
In this study, stiffness estimates were calculated for normal connective tissues and for a condition in which intermuscular connectivity was enhanced. Such a condition simulates the effects of altered connective tissues as a consequence of trauma, invasive surgical interventions, muscle-tendon injury, aging and pathologies such as muscle spasticity affecting connective tissues. For instance, scar tissue formation following rupture of muscle fibers and of the surrounding connective tissues framework has been associated with shearing type of muscle injuries (Kääriäinen et al. 2000). Scar tissue may result also in a reduction of the over-all muscle belly displacement and a localized increase of tissue strain (Silder et al. 2010). We have shown that such enhancement of intermuscular connections alters the mechanical interaction between synergistic muscles (Bernabei et al. 2016). Similar mechanical effects of scar tissue may be expected following orthopedic surgical interventions involving disruption of muscle fascia or other connective tissues (e.g. fasciotomy, aponeurotomy and tenotomy). In fact, several studies have reported changes of intermuscular interaction following resection of connective tissues that served as pathways of force transmission (Smeulders et al. 2002;Maas et al. 2013;Ateş et al. 2013). A first attempt to model effects of altered non-myotendinous force transmission due to scar tissue development after tendontransfer surgery was reported by Andersen (2009). In his extensive musculoskeletal model, a constant factor describes how muscle force is split between the transferred rectus femoris muscle (RF) and the vasti muscles via scar tissue. However, fibrotic connective tissues were defined as spurious muscle elements capable of producing a percentage of the normal muscle strength, regardless of the relative position between RF and patella. In this respect, estimates of inter-and extramuscular stiffness could complement Andersen's model by taking into account a more realistic nonlinear stiffness factor, dependent on muscles' relative position and scaled to the extent of tissue-integration post-intervention (see Fig. 5b). As opposite to the increased force transmission with stiffer epimuscular tissues, a degradation of mechanical properties of connective tissues has been proposed to explain reduced muscle force in Ehlers-Danlos patients (Huijing et al. 2010). Effects of degraded mechanical properties of connective tissues, e.g. a lower stiffness, could be extrapolated from the model presented here, which would provide a qualitative description of the tendon force changes as a consequence of impaired force transmission. The complex organization of the connective tissues network surrounding muscles requires a detailed description of the anatomy of force pathways. Recently, subject-specific models to predict effects of orthopedic surgery have been proposed to overcome variability among patients in terms of anatomy, activity levels, loading conditions etc. (Carbone et al. 2015). In this regard, the approach described here would help obtaining more accurate estimates of human function when connective tissues pathways are affected or impaired after muscle-tendon injury or orthopedic surgery. Further development is required to scale the non-myotendinous stiffness estimates to human muscles, also taking into account the organization and heterogeneity of inter-and extramuscular linkages between human muscle groups. In addition, the magnitude of non-myotendinous force transmission appears to be related to the extent of muscle activation and coactivation as suggested by differences between passive and active states (Finni et al. 2015;Tijs 2015). Therefore, scaling of force transmission estimates with muscle activation levels may be needed to integrate this approach into a musculoskeletal model driven by neuromuscular activation time series.

Conclusions
The multi-muscle phenomenological model presented here provides a tool for studying effects of epimuscular force transmission on synergistic muscles force production. It improves our representation of muscle mechanics by implementing the connective tissues network, which surrounds muscles and can transmit the force produced within the muscle to other muscles and to non-muscular structures. Thus, this approach enables predictions of tendon forces in musculoskeletal models when the common assumption of muscle mechanical independence is challenged. Direct applications involve simulation of scar tissue development and prediction of surgical outcomes, where the effects of altered epimuscular myofascial force transmission may be substantial.
The stiffness estimate assessed in this study is a lumped representation of the force redirected from one muscle to the tendon of a neighboring muscle via intermuscular connective tissues (INT), blood vessels and nerves (NV1, NV2). The constituent equation for the intermuscular and extramuscular connectivity elements was defined as: where F is force change, P R is the change in relative position between connected muscles andK is the parameter describing the stiffness of a given force pathway. Connective tissues, ranging from tendon to skin to blood vessels, tend to have quasi-static mechanical properties: the stiffness increases fairly linearly with force over the primary operating range (Fung 1981;Winters and Woo 1990). Thus, the lumped stiffness of intermuscular connections was modeled as dependent on the synergists' relative position P R for which a piecewise linear relation was chosen (Eq. 7). Each muscle was divided into a proximal and a distal node to conveniently represent myofascial connective tissue anchoring points (Fig. 1c).

Model assumptions
(a) Muscle bellies were assumed to undergo negligible deformation during contraction. (b) Sarcomere length changes were assumed to be negligible as muscles were kept at constant MTU length, so that constant isometric force was obtained at each P R . It should be noted that sarcomere length distribution may vary, at least in passive conditions, for different imposed relative positions between muscles (Tijs et al. 2015b). (c) Length changes of underlying structures (INT, NV1, NV2 pathways) were considered proportional to muscles' relative displacement: (10) (f) The ratio between length changes of NV1 before (intact) and after resection of the myofascial linkages (post I) was assumed to be a constant: Sensitivity analysis on c was performed by calculating the residuals between the range of predicted force transmission values F INT , as defined in Eq. 21, for c ∈ (0.05; 1.00) and the F INT obtained by setting c = 0.9, corresponding to minimal changes of l NV1 after resection of INT (post I). Analysis of root-mean-square of residuals revealed that the difference obtained by assuming c = 0.9 was up to 6.0 % of the average F INT over the imposed range P R = 6mm.

Stiffness estimates assessment
With intact intermuscular and extramuscular connective tissue pathways, the balance of forces produced by active muscles within the system shown in Fig. 1 can be written as: The aforementioned conditions assessed in the in situ experiment were: (i) intact connective tissues (intact), (ii) full resection of intermuscular myofascial connective tissues (INT resection, post I) and (iii) removal of all connections between SO and LG + PL (INT and NV1 resection, post II). Given the isometric conditions, accelerations were inherently equal to zero (u LG P L prox =u LGPL prox ), therefore: After resection of the intermuscular connective tissues (the INT pathway), this becomes The residual LG + PL proximal-distal force difference (F post I I NV2 ) with P R ∈ (−3; +3) mm was measured and interpolated with a III-order polynomial (Fig. 7, dotted line). Considering Eqs. 9, 15 and Eq. 16 becomes: By adding this terms to the force equilibrium before resection (intact condition), Eq. 14 becomes:  LG + PL proximal-distal force difference ( F LGPL ) with progressive resection of connective tissues pathways. F LGPL with intact connective tissues (intact), after resection of intermuscular connective tissues (post I) and after resection of all intermuscular and extramuscular force pathways (post II), but the neurovascular tract running proximally from LG + PL (NV2). Residual LG + PL proximal-distal force difference data points over the applied range of SO/LG + PL relative displacement after resection of intermuscular pathways and the neurovascular tract running between SO and LG + PL (post II resection) were interpolated with a III-order polynomial (dashed line, n = 1) which considering Eq. 11, can be written as: LGPL + F post II LGPL · (c − 1) (20) The force transmitted via intermuscular connective tissue pathways (INT) between SO and LG + PL as a function of synergists' relative position (P R ) can now be written as LGPL − c * F post I LGPL + F post II LGPL * (c − 1) (21) Thus, according to Eq. 7, the estimate of stiffness associated with intermuscular myofascial linkages was defined aŝ While the force transmitted via the intermuscular neurovascular tract between SO and LG+PL (NV1) as a function of synergists' relative position change (P r ) can be written as F NV1 (P r ) = F post I LGPL − F post II LGPL Fig. 8 Effects of distal versus full resection on the estimate of intermuscular force transmission. Intermuscular force estimates (mean ± SD) over a range of muscle relative displacement obtained by comparing LG+PL tendon forces before and after resection of all intermuscular connective tissues (full resection, white circles) or distal intermuscular connective tissues only (distal resection, black circles). The ratio between full resection and distal resection force values provided the r scaling factor at each considered P R Thus, the relative estimate of stiffness was defined aŝ

Model validation
The validation of the model obtained by the calibration dataset was performed by comparing predictions of F LGPL ( F mod LGPL ) with experimentally measured values of F LGPL ( F exp LGPL ), which belong to a testing dataset. Such dataset was obtained from the same animals, but applying a different protocol that spans physiological muscle lengths and relative positions changes. The predicted LG + PL force difference which describes the amount of force redirected inter-and extramuscularly was defined as: We expected that, for physiological muscle lengths and relative position of SO/LG+PL, a more limited involvement of intermuscular and extramuscular pathways would result in a reduced magnitude of force transmission (Bernabei et al. 2015). This was taken into account by scaling the stiffness descriptor of INT pathway with a function of the proximal and distal relative displacement between muscles: r ( L R ). This function was determined experimentally by resecting only the distal portion of the intermuscular connective tis-sues and calculating the ratio of changes in force between full resection and distal resection for each relative position (Fig. 8, n = 3). The implicit assumption is that the absence of distal connections would resemble the condition in which only the proximal connections are involved in the force transmission pathway. Therefore, with a proximal lengthening of LG + PL muscle that corresponds to physiological conditions, the predicted LG + PL force difference was calculated as: