1 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 extra- and 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. 2000; Smeulders 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 (\(\hbox {LG}+\hbox {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 \(\hbox {LG}+\hbox {PL}\) for both normal and enhanced intermuscular connectivity. Finally, the relative contribution of extramuscular pathways of force transmission, i.e., SO and \(\hbox {LG}+\hbox {PL}\) neurovascular tracts, with respect to the epimuscular myofascial and myotendinous pathways was assessed.

2 Methods

2.1 In situ experiment

2.1.1 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 \(\hbox {LG}+\hbox {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.

2.1.2 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 \(\hbox {LG}+\hbox {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 \(\hbox {LG}+\hbox {PL}\), a tissue-integrating mesh (\(n=7\); Premilene\(^{\circledR }\) mesh, B. Braun Melsungen AG, Germany) was sutured to the SO muscle belly at the interface of SO and \(\hbox {LG}+\hbox {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\(^{\circledR }\), 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.

2.1.3 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 \(\hbox {LG}+\hbox {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 \(\hbox {LG}+\hbox {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.

Fig. 1
figure 1

Schematic representation of the experimental setup and of the model. a Epimuscular (inter- and extramuscular) pathways which were left intact after dissecting the plantar flexors synergistic group before the in situ measurement. A proximal and a distal intermuscular pathway of force transmission were distinguished, while the connective tissue bundles enwrapping blood vessels and nerves proximal to \(\hbox {LG}+\hbox {PL}\) (NV2) and at the \(\hbox {SO}/\hbox {LG}+\hbox {PL}\) interface (NV1) were defined as extramuscular pathways of force transmission. b Proximal and distal relative displacement \((P_{R})\) between SO and \(\hbox {LG}+\hbox {PL}\) MTUs muscles imposed during the in situ experiment. Tendon forces were measured at the proximal and distal \(\hbox {LG}+\hbox {PL}\) tendons as well as at the distal SO tendon. The proximal SO tendon insertion on the tibia was left intact. c Lumped-parameter model rendering SO and \(\hbox {LG}+\hbox {PL}\) individual actuators with a simple Hill-model composed by a contractile element (CE), both series (SE) and parallel (PE) elastic elements and proximal and distal mass-nodes \((\hbox {m}_{\mathrm{PROX, DIST}})\). The mechanical interaction provided by intermuscular \((\hbox {INT}_{\mathrm{PROX, DIST}})\) and extramuscular (NV1, NV2) pathways of force transmission were represented with piecewise linear elastic elements between SO and \(\hbox {LG}+\hbox {PL}\) lumped-mass elements. \(K_{\mathrm{SO}}^{T}, K_{\mathrm{LGPT}}^{T}\): tendon stiffness. \(K_\mathrm{NV1}, K_\mathrm{NV2}\): Neuro-vascular tract stiffness (extramuscular), \(K_{\mathrm{INT}}\); intermuscular myofascial connective tissues stiffness. \(I, \Delta l\): length, length changes in the given element. CE, PE, SE: contractile element, parallel elastic element series elastic element. \(\Delta F_{\mathrm{LGPL}}\): proximal-distal \(\hbox {LG}+\hbox {PL}\) force difference (intact: intact condition; post I: myofascial pathway severed; post II: intermuscular neurovascular tract severed). \(P_{R}\): positions of \(\hbox {LG}+\hbox {PL}\) relative to SO muscle, \(P_{\mathrm{REF}}\): \(P_{R}=0\,\hbox {mm}\). Muscles: reference position corresponding to \(90^{\circ }\)\(90^{\circ }\) knee-ankle joint angles. \(\widehat{{F}_\mathrm{INT}},\widehat{{F}_\mathrm{NV1}}, \widehat{{F}_\mathrm{NV2}}\) :force estimates of non-myotendinous pathways of force transmission, calculated at each different SO/LG+PL relative position \(P_{R}\) from experimental data. \(\hat{K}_{\mathrm{INT}}\), \(\hat{K}_{\mathrm{NV1}}\), \(\hat{K}_{\mathrm{NV2}}\): stiffness estimates of non-myotendinous pathways of force transmission; intermuscular \(({ INT})\) and extramuscular (\({ NV1}, { NV2}\))

2.1.4 Nerve stimulation

SO and \(\hbox {LG}+\hbox {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.

2.1.5 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^{\circ }\), which was defined as the reference position \((P_{\mathrm{REF}})\). Markers were placed near the SO distal and \(\hbox {LG}+\hbox {PL}\) proximal and distal tendons to apply different relative positions between SO and \(\hbox {LG}+\hbox {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 \(\hbox {LG}+\hbox {PL}\) relative to SO muscle \((P_{R})\). The proximal and distal \(\hbox {LG}+\hbox {PL}\) tendons were repositioned in steps of 1 mm from \(\hbox {P}_{\mathrm{REF}}-3\,\hbox {mm}\) in distal direction to \(\hbox {P}_{\mathrm{REF}}+3\,\hbox {mm}\) in proximal direction (Fig 1b). \(\hbox {LG}+\hbox {PL}\) tendons were repositioned together, while the distal and proximal tendons of SO were kept at the same position corresponding to \(P_{\mathrm{REF}}\). Thus, MTU lengths were constant over all the imposed conditions, with SO and \(\hbox {LG}+\hbox {PL}\) operating on the ascending region of their length–force curves (mean ± SD at \(\hbox {P}_{\mathrm{REF}}: \hbox {F}_{\mathrm{LG+PL}}=11.5\pm 0.2\hbox {N}, \hbox {F}_{\mathrm{SO}}= 1.3\pm 0.5\hbox {N}\), 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 \(\hbox {LG}+\hbox {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).

Fig. 2
figure 2

Exemplar waveforms of forces exerted at the tendons of SO and \(\hbox {LG}+\hbox {PL}\). Passive and total forces at soleus distal (SO), \(\hbox {LG}+\hbox {PL}\) proximal (\(\hbox {LG}+\hbox {PL}\) prox) and \(\hbox {LG}+\hbox {PL}\) distal (\(\hbox {LG}+\hbox {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 \(\hbox {LG}+\hbox {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 \(\hbox {LG}+\hbox {PL}\) muscle fibers

2.1.6 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 \(\hbox {LG}+\hbox {PL}\) tendons simultaneously, we calculated (i) changes in force exerted at the distal tendon of SO \((\Delta F_\mathrm{SO})\), with SO kept at a constant length and (ii) changes in the force difference between the proximal and distal tendons of \(\hbox {LG}+\hbox {PL}\), with \(\hbox {LG}+\hbox {PL}\) kept at a constant length. Changes in SO force were expressed relative to the value at \(P_\mathrm{REF} (\Delta F_\mathrm{SO}=F_\mathrm{SO} (P_{R})-F_\mathrm{SO} (P_\mathrm{REF}))\). The \(\hbox {LG}+\hbox {PL}\) force difference \((\Delta F_\mathrm{LGPL})\) was calculated by subtracting the force exerted distally from the force exerted proximally (i.e., \(F_\mathrm{PROXIMAL}-F_\mathrm{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 (Huijing and Baan 2001; Maas and Sandercock 2010). \(\Delta F_\mathrm{LGPL}\) and \(\Delta F_\mathrm{SO}\) as defined above will be referred to as the calibration dataset.

To test for effects of muscle relative position on proximal and distal \(\hbox {LG}+\hbox {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 \(\Delta F_\mathrm{SO}\) and \(\Delta F_\mathrm{LGPL}\) over \(P_{R} \in [-3;+3]\,\hbox {mm}\). \(\Delta F_\mathrm{SO}\) and \(\Delta F_\mathrm{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.

2.2 Modeling study

2.2.1 Inter- and extramuscular pathways of force transmission

The myofascial connective tissues at the interface between SO and \(\hbox {LG}+\hbox {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.

2.2.2 Model description

A phenomenological lumped-parameter model was used to represent SO and \(\hbox {LG}+\hbox {PL}\) interactions (Fig. 1c). Each muscle 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 inter- and extramuscularly to surrounding structures (Maas et al. 2004). Thus, \(\Delta F_\mathrm{LGPL}\) was considered as an estimate of the force redirected from \(\hbox {LG}+\hbox {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_\mathrm{proximal}\) and \(F_\mathrm{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:

$$\begin{aligned}&\hbox {(node 1)} \qquad F_\mathrm{LGPL}^\mathrm{prox} -F_\mathrm{NV2}-F_\mathrm{LGPL}^\mathrm{PE}-F_\mathrm{LGPL}^\mathrm{CE} =0 \end{aligned}$$
(1)
$$\begin{aligned}&\hbox {(node 2)}\qquad F_\mathrm{LGPL}^\mathrm{dist}-F_\mathrm{LGPL}^\mathrm{PE}-F_\mathrm{LGPL}^\mathrm{CE} +F_\mathrm{INT} + F_\mathrm{NV1}=0 \end{aligned}$$
(2)

where F is the force redirected to the given node via the superscripted element. Because MTUs of \(\hbox {LG}+\hbox {PL}\) and SO were kept at the same length between conditions and muscle contractions were isometric, \(F_\mathrm{CE}\) and \(F_\mathrm{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_\mathrm{INT}\). Solving Eqs. 12 for \(F_\mathrm{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).

2.2.3 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 \((\hat{K})\) for each inter- and extramuscular pathway of force transmission (INT, NV1, NV2) was defined as the ratio of force increase \((\Delta F)\) over a unit of relative displacement \((\Delta 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 \(\hbox {SO}/\hbox {LG}+\hbox {PL}\) relative displacement were calculated as:

$$\begin{aligned} \widehat{F_{\mathrm{INT}}}(P_{r})= & {} \Delta F_\mathrm{LGPL}^{\mathrm{intact}}-c*\Delta F_\mathrm{LGPL}^{\mathrm{post}\,\mathrm{I}} +\Delta F_\mathrm{LGPL}^{\mathrm{post}\,\mathrm{II}}*\left( {c-1} \right) \qquad \mathrm{[N]} \end{aligned}$$
(3)
$$\begin{aligned} \widehat{F_\mathrm{NV1}}(P_{r})= & {} \Delta F_\mathrm{LGPL}^{\mathrm{post}\,\mathrm{I}} -\Delta F_\mathrm{LGPL}^{\mathrm{post}\,\mathrm{II}} \qquad \hbox {[N]} \end{aligned}$$
(4)
$$\begin{aligned} \widehat{F_\mathrm{NV2}}(P_{r})= & {} \Delta F_\mathrm{LGPL}^{\mathrm{post}\,\mathrm{II}}\qquad \hbox {[N]} \end{aligned}$$
(5)

where \(\Delta F_\mathrm{LGPL}\) is the force difference at \(\hbox {LG}+\hbox {PL}\) proximal-distal 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\in (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 \(\hbox {LG}+\hbox {PL}\) proximal tendon force were calculated \((\widehat{F_{\mathrm{INT}}}/F_\mathrm{LGPL}^\mathrm{prox}\,[\%]\) and \(\widehat{F_\mathrm{NV1}}/F_\mathrm{LGPL}^\mathrm{prox}\, [\%], \widehat{F_\mathrm{NV2}}/F_\mathrm{LGPL}^\mathrm{prox}\, [\%]\)).

2.2.4 Model validation

The lumped-parameter model was validated by comparing the predicted \(\hbox {LG}+\hbox {PL}\) tendon force balance \((\Delta F_\mathrm{LGPL}^{\mathrm{mod}})\) with previously collected experimental data (Bernabei et al. 2016), here referred as the testing dataset \((\Delta F_\mathrm{LGPL}^\mathrm{exp})\). The experimental protocol used for measuring forces in the testing dataset differed from that used for the calibration dataset. For the testing dataset, \(\hbox {LG}+\hbox {PL}\) force changes were measured applying a physiological range of lengths and relative displacements \((P_{R}^\mathrm{prox})\) between SO and \(\hbox {LG}+\hbox {PL}\), with only \(\hbox {LG}+\hbox {PL}\) proximal tendon repositioned from \(P_{R}^\mathrm{prox}-3\,\hbox {mm}\) to \(P_{R}^\mathrm{prox}+3\,\hbox {mm}\), while the distal tendons were kept at the same position. For the calibration dataset instead, both proximal and distal tendons of \(\hbox {LG}+\hbox {PL}\) were repositioned together from \(P_{R} - 3\,\hbox {mm}\) to \(P_{R} +3\,\hbox {mm}\), with constant muscle lengths for \(\hbox {LG}+\hbox {PL}\) and SO (see above). Predicted values of \(\hbox {LG}+\hbox {PL}\) proximal-distal force difference were calculated as follows:

$$\begin{aligned} \Delta F_\mathrm{LGPL}^\mathrm{mod}(\Delta P_{R}^\mathrm{prox})= & {} \Delta F_\mathrm{LGPL}^\mathrm{mod}(P_\mathrm{REF})+r\cdot {\hat{K}}_\mathrm{INT}\Delta P_R^\mathrm{prox}\nonumber \\&+\,\hat{K}_\mathrm{NV1}\cdot \Delta P_{R}^\mathrm{prox}\nonumber \\&+\,\hat{K}_\mathrm{NV2}\cdot \Delta P_R^\mathrm{prox} \qquad \hbox {[N]} \end{aligned}$$
(6)

where \(P_\mathrm{REF}\) corresponded to \(P_{R}=0\); this \(\hbox {SO}/\hbox {LG}+\hbox {PL}\) relative position represented the matching condition for muscle lengths and relative displacement between the calibration dataset and the testing dataset. \(\Delta P_{R}^\mathrm{prox} \) defined the proximal displacement of \(\hbox {LG}+\hbox {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.

2.2.5 Analysis of model results

ANOVA for repeated measurements (within factor: \(\hbox {LG}+\hbox {PL}\) proximal tendon position, between factor: predicted vs. measured data) was used to test for differences between the values of \(\hbox {LG}+\hbox {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 \(\hbox {LG}+\hbox {PL}\) force differences. The bias was calculated as the mean difference between measured \((\Delta F_\mathrm{LGPL}^\mathrm{exp})\) and predicted force values \((\Delta F_\mathrm{LGPL}^\mathrm{mod})\).

3 Results

3.1 Regression analysis

3.1.1 Tetanic forces

A linear relationship between \(\Delta F_\mathrm{LGPL}\) and \(\Delta F_\mathrm{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 \(\hbox {LG}+\hbox {PL}\) total forces suggests that the force discrepancy at \(\hbox {LG}+\hbox {PL}\) tendons was balanced by the force exerted at SO distal tendon, which can be explained by force transmitted from \(\hbox {LG}+\hbox {PL}\) to SO via connective tissues at their muscle belly interface.

Fig. 3
figure 3

Comparison of SO and \(\hbox {LG}+\hbox {PL}\) force discrepancies in active and passive state of muscles. Regression analysis between \(\Delta F_{\mathrm{SO}}\) (mean \(+\) SD) and \(\Delta F_{\mathrm{LGPL}}\) (mean \(+\) SD) with \(\hbox {SO}/\hbox {LG}+\hbox {PL}\) relative displacement from \(P_{R}=-3\) to \(P_{R}=+3\,\hbox {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

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: \(\Delta F_\mathrm{LGPL}=0.51\hbox {N}\), TI: \(\Delta F_\mathrm{LGPL}=2.12\hbox {N}\) at \(P_{R}=+3\,\hbox {mm}\)). Finally, it should be noted that \(\Delta F_\mathrm{LGPL}\) reached 0 N at \(P_{R}=0\hbox {mm}\) (REF position) and became negative for \(P_{R}\in \) \([-3; 0]\) mm, which suggests a reversal of the sign of the epimuscular force transmission vector at \(P_{R}=0\;\hbox {mm}\).

Fig. 4
figure 4

Epimuscular force transmission with normal and enhanced connective tissues. Distribution of epimuscularly transmitted force among intermuscular connective tissues (INT) and extramuscular neurovascular tracts (NV1, NV2), over each step of \(\hbox {SO}/\hbox {LG}+\hbox {PL}\) MTUs displacement (\(\Delta P_{R}\)). Force values are expressed as a percentage relative to the \(\hbox {LG}+\hbox {PL}\) proximal tendon total force (\(\widehat{F}/F_{\mathrm{prox}}^{\mathrm{LG}}\) [%]) 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 \(\hbox {LG}+\hbox {PL}\) tendon, is shown for \(P_{R}=2\,\hbox {mm}\)

3.1.2 Passive forces

The relationship between passive \(\Delta F_\mathrm{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 \(\hbox {LG}+\hbox {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 \(\Delta F_\mathrm{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.

3.2 Estimates of force transmission via intermuscular and extramuscular pathways

Force generated by \(\hbox {LG}+\hbox {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 \((\hat{F}/F_\mathrm{LGPL}^\mathrm{prox})\) for NO (Fig. 4a) and TI (Fig. 4b) was maximal at \(P_{R}=-3\,\hbox {mm}\) (7.7 %) and at \(P_{R}=+3\,\hbox {mm}\) (17.7 %), respectively. NV1 and NV2 contributed to transmit force up to 1.1 and −4.0 % of proximal \(\hbox {F}_{\mathrm{LGPL}}\) at \(P_{R}=-3\,\hbox {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 \(\hbox {LG}+\hbox {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 configuration, i.e., angle and direction of the force transmission vector in response to changes in SO/\(\hbox {LG}+\hbox {PL}\) relative position.

3.3 Validation study

ANOVA showed no significant differences between \(\Delta F_\mathrm{LGPL}^\mathrm{exp}\) and \(\Delta F_\mathrm{LGPL}^\mathrm{mod}\) (\(p=.956\)). In the testing dataset, \(\hbox {LG}+\hbox {PL}\) proximal-distal force difference varied with \(\hbox {LG}+\hbox {PL}\) proximal position \(\Delta P_R^\mathrm{prox} \) 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 \(\Delta F_\mathrm{LGPL}\) for NO and underestimated at high \(\Delta F_\mathrm{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).

Table 1 Stiffness estimates (calibration dataset)

4 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 muscles. 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 (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).

Fig. 5
figure 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 \(\hbox {LG}+\hbox {PL}\) were applied (Bernabei et al. 2016). Averaged predicted values (black circles) of \(\hbox {LG}+\hbox {PL}\) proximal-distal force difference with repositioning \(\hbox {LG}+\hbox {PL}\) proximal tendon \((\Delta P_R^{\mathrm{prox}})\) from \(-3\) to \(+3\) mm from \(P_{\mathrm{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 connectivity (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 \(\Delta F_{\mathrm{LG}+\mathrm{PL}}\) on the x-axis were expressed as a percentage of the total range of \(\hbox {LG}+\hbox {PL}\) proximal-distal force difference \((\max (\Delta F_{\mathrm{LG}+\mathrm{PL}})-\min (\Delta F_{\mathrm{LG}+\mathrm{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)

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 (Huijing and Baan 2001; 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 (Huijing and Baan 2003; 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 \(\hbox {LG}+\hbox {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 force transmission phenomena, which could be integrated in a more comprehensive musculoskeletal model. This would enable predictions of force changes at the tendons of non-isolated muscles in functional tasks when connective tissues are altered.

Fig. 6
figure 6

Epimuscular force transmission with knee angle extension. Force transmission estimates predicted by the model \((\Delta F_{\mathrm{LGPL}}^{\mathrm{mod}})\) plotted a as a function of knee angle corresponding to the applied proximal \(\hbox {LG}+\hbox {PL}\) length changes \((P_R^{\mathrm{prox}})\) 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^{\circ }\) knee angle are highlighted for the evaluated intermuscular connectivity conditions, i.e., normal (NO) and enhanced connectivity (TI). The secondary axis expresses the predicted \(\Delta F_{\mathrm{LGPL}}^{\mathrm{mod}}\) as a percentage of the mean \(\hbox {LG}+\hbox {PL}\) proximal tendon force (\(\hbox {F}_{\mathrm{LG}+\mathrm{PL}}\), proximal: \(11.5\pm 0.2\hbox {N}\)). In the bottom panel, linear interpolation of force transmission estimates between NO (green) and TI (red) over the range of corresponding knee angles

4.1 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 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 \(\hbox {LG}+\hbox {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 \(\Delta F_\mathrm{LGPL}\) and \(\Delta F_\mathrm{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.

4.2 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 large-scale 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.

4.3 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 overall 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 tendon-transfer 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 co-activation 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.

5 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.