Modelling motor units in 3D: influence on muscle contraction and joint force via a proof of concept simulation

Functional heterogeneity is a skeletal muscle’s ability to generate diverse force vectors through localised motor unit (MU) recruitment. Existing 3D macroscopic continuum-mechanical finite element (FE) muscle models neglect MU anatomy and recruit muscle volume simultaneously, making them unsuitable for studying functional heterogeneity. Here, we develop a method to incorporate MU anatomy and information in 3D models. Virtual fibres in the muscle are grouped into MUs via a novel “virtual innervation” technique, which can control the units’ size, shape, position, and overlap. The discrete MU anatomy is then mapped to the FE mesh via statistical averaging, resulting in a volumetric MU distribution. Mesh dependency is investigated using a 2D idealised model and revealed that the amount of MU overlap is inversely proportional to mesh dependency. Simultaneous recruitment of a MU’s volume implies that action potentials (AP) propagate instantaneously. A 3D idealised model is used to verify this assumption, revealing that neglecting AP propagation results in a slightly less-steady force, advanced in time by approximately 20 ms, at the tendons. Lastly, the method is applied to a 3D, anatomically realistic model of the masticatory system to demonstrate the functional heterogeneity of masseter muscles in producing bite force. We found that the MU anatomy significantly affected bite force direction compared to bite force magnitude. MU position was much more efficacious in bringing about bite force changes than MU overlap. These results highlight the relevance of MU anatomy to muscle function and joint force, particularly for muscles with complex neuromuscular architecture. Supplementary Information The online version contains supplementary material available at 10.1007/s10237-022-01666-2.


Introduction
The motor unit is the fundamental element of force production in the neuromuscular system, comprising all muscle fibres that are innervated by a single motor neuron (Sherrington 1929). A motor unit is typically classified by the number of fibres that belong to it (size), and the contractile properties of those fibres. For example, type-S, -FR, and -FF motor units comprise, type I, type IIa, and IIb fibres, respectively. The number of recruited motor units together with their rate of stimulation determines the force exerted by a muscle (Adrian and Bronk 1929). In general, motor unit recruitment follows Henneman's size principle (Henneman 1957), although it can be altered by the type of task, contraction velocity, or sensory feedback (Schindler et al. 2014;Ogawa et al. 2006;Romaiguére et al. 1989;Hodson-Tole and Wakeling 2009).
Besides varying in size, motor units also show different anatomy. Motor unit anatomy refers to the distribution of its fibres across the muscle and the cross-sectional area enclosing these fibres refers to its territory. While muscle anatomy varies greatly both within and between individuals, motor unit territories show common features such as size variation, overlap with other territories, and irregular shape (Bodine-Fowler et al. 1990;. Additionally, 1 3 motor unit types may be restricted to certain parts of the muscle cross-sectional area (Mesin et al. 2010;Tonndorf and Hannam 1994). The deep portions of the tibialis anterior, for example, have a higher proportion of type-FF and -FR motor units (Henriksson-Larsén et al. 1985;Mesin et al. 2010).
The interplay between the selective recruitment of motor units and their anatomy can lead to heterogeneous muscle contraction. This phenomenon was described by Herring et al. as "a single muscle's capacity to create a multitude of force vectors through selective localised contraction" (Herring et al. 1979). The temporalis and masseter motor unit, for example, are selectively recruited depending on bite force direction (Blanksma et al. 1990;Ogawa et al. 2006). (Blanksma et al. 1990) used fine-wire electrodes to measure activity in different regions of the temporalis and found that, apart from the anterior regions, the muscle showed alterations in activity in response to changes in bite force direction. (Phanachet et al. 2003) used fine-wire electrodes to record electromyography signals of single motor unit at various locations in lateral pterygoid and found that 14 % of motor units were specialised to one task. Similarly, (Schindler et al. 2014) used fine-wire electrodes in the masseter and found that about half (46 %) of the investigated motor unit were specialised to one task. (Staudenmann et al. 2009) used a grid of surface electromyography electrodes to record the activity of the triceps surae during plantar flexions along predefined directions while recording forces at the foot and showed correlations between different surface activity and the directions of recorded force. Investigating the same muscle (group), (Wakeling 2009) used a coarser surface electromyography grid to measure muscle activity during cycling against different resistances and showed very little correlation within regions of the same muscle, implying selective recruitment of the muscle regions. These studies, among others (ter et al. 1984;Pérot et al. 1996;Turkawski et al. 1998;Guzmán et al. 2015;Staudenmann et al. 2009;Murray et al. 1999;Blanksma et al. 1955), suggest that the task dependence of skeletal muscles is the norm rather than the exception.
Human movement is often investigated by measuring the body's kinematics and estimating corresponding muscle activations by applying inverse-dynamics to biomechanical models (e.g. Lloyd and Besier 2003;Pizzolato et al. 2017;Harandi et al. 2020). Typically, such models are based on multi-body simulations, which represent muscles as lumped contractile units, i.e. they neglect the functional organisation of skeletal muscles into motor units. At most, only statistical information about motor unit firing is used (e.g. Sartori et al. 2017). Therefore, the interplay between motor unit anatomy and human movement remains largely unexplored. In contrast, 3D, macroscopic, continuum-mechanical muscle models (3D models) do indeed have the ability to spatially resolve motor units. However, such models have largely ignored this feature of muscle contraction (e.g. Johansson et al. 2000;Fernandez et al. 2005;Blemker and Delp 2005;Röhrle and Pullan 2007;Wu et al. 2014;Weickenmeier et al. 2017;Chi et al. 2010;Péan et al. 2019)-most likely due to the representation of muscle fibres as a continuous vector field, which makes the grouping of fibres to define motor unit anatomy challenging. Currently, only multi-scale continuum-mechanical muscle models resolve motor unit anatomy by describing individual muscle fibres (e.g. Heidlauf and Röhrle 2014;Schmid et al. 2019;Röhrle et al. 2012;Marcucci et al. 2017;Teklemariam et al. 2019). However, since these multi-scale models are computationally very expensive to solve, existing simulation studies have been restricted to idealised, isolated muscles or fibre bundles.
Within this work, we propose an approach that combines both motor unit activity and anatomy into 3D models to represent the heterogeneous nature of muscle contraction. We do this by decomposing muscle activity into motor unit activity and anatomy, resulting in a spatiotemporal muscle activation parameter. This requires the computation of physiologically plausible motor unit distributions, activity, and a method to incorporate this information into a macroscopic continuum-mechanical framework. First, a brief overview of continuum-mechanical theory and constitutive modelling is provided (Sect. 2.1). Followed by the methods for generating motor unit information. Motor unit anatomy is created via a novel "virtual innervation" technique (Sect. 2.2.1), while motor unit activity models are drawn from the literature (Sect. 2.2.2). Discrete motor unit anatomy, on the muscle fibre scale, is then incorporated into the continuum mechanical, finite element framework, requiring a mapping of the (microscale) motor unit properties to the (macroscopic) finite element mesh and an assumption on action potential conduction velocity. An idealised 2D geometry is used to investigate the mesh dependency of muscle activity computation (Sect. 2.3.1) and an idealised 3D geometry to verify the influence of instantaneous action potential propagation on force production (Sect. 2.3.2). Finally, we demonstrate virtual innervation on the masseters of a masticatory model and generate a range of unique motor unit distributions. Each unique masseter is activated by an identical motor unit recruitment protocol to simulate a static bite. The differences in the bite force reveal the interplay between motor unit distribution and muscle function (Sect. 2.3.3). Our technique for creating motor unit-enriched 3D models provides a transparent environment for studying relationships between motor unit physiology, muscle function, and movement, such as muscle task dependency.

Three-dimensional musculotendon modelling
Skeletal muscles undergo large deformations during contraction and movement, making the theory of finite elasticity a suitable choice to describe their mechanical response. Here, we provide a very brief sketch of the theory. For further details, see Holzapfel (2000); Bonet and Wood (1997) for example.
Consider the musculotendon complex as a continuum body B 0 (at t = 0 , B otherwise) comprised of material points P . Each material point P has coordinates X in the reference configuration ( t = 0 ), being mapped to x in the current configuration (at time t). This mapping is described by the placement function x = (X, t) . Then, the deformation and motion of each material point P is described by the deformation gradient tensor F = x∕ X , and the right Cauchy-Green deformation tensor C = F T F provides a rotation invariant deformation measure.
The balance of linear momentum governs the deformation of B and is given here in its local form with respect to the current configuration, i.e.
where u = x − X is the displacement vector in m, is the mass density in kg/m 3 and is the Cauchy stress tensor in N/m 2 . Body forces, such as gravity, are absent in the momentum balance as we assume them to be negligible compared to muscle forces.
Constitutive modelling The stress in Eq. 1 describes the material's mechanical behaviour based on its constitutive relation. Details of the particular constitutive relation used to describe the musculotendon complex are given in Saini and Röhrle (2022);  and are only briefly covered here. We assume that skeletal muscle is hyperelastic, transversely isotropic, and nearly incompressible. Furthermore, the passive and active stresses are additively decomposed, i.e. S = S pass + S act , where S is the second Piola-Kirchhoff stress tensor. The passive stress is further decomposed between contributions from the (isotropic) bulk and along the (anisotropic) fibre direction, i.e. S pass = S iso + S aniso .
The passive stresses are derived from strain energies, with S iso and S aniso characterised by Mooney-Rivlin (Mooney 1940;Rivlin and Rideal 1948) and Balzani et al. (206) formulations, respectively. Accordingly, the isotropic stress is described by and the anisotropic stress by (1) div =ü, In Eqs. 2 and 3, c 1−6 are independent material parameters, all in N/m 2 except for c 4 and c 6 , which are dimensionless. The parameter d = 2(c 1 + 2 c 2 ) ensures a stress-free reference configuration, and c penalises dilatational (volume changing) deformations. The Jacobian J = det(F) deviates from unity under such dilatational deformations. The invariants of C are defined by I 1 = tr(C) and I 4 = tr(MC) = 2 f , where f is the fibre stretch. Lastly, M = a 0 ⊗ a 0 is a structural tensor, which describes the (transverse) anisotropy of the material via the fibre direction in the reference configuration a 0 .
The active stress is described by where S max is the peak (isometric) stress that the muscle can produce, in N/m 2 , and f l ∈ [0, 1] represents the normalised amount of (actin and myosin) filament overlap, i.e. the forcelength relationship of muscle. It is characterised by where opt f is the optimal fibre length at which S max is produced. The parameters i and v i determine the shape of the force-length curve, where the subscripts i = asc and i = dsc when f < opt f and f > opt f , representing the ascending and descending branches of the force-length relationship, respectively.
Lastly, m ∈ [0, 1] represents the normalised muscle activity, varying between 0 and 1 for fully passive and active states, respectively. Traditionally this is only a temporally varying quantity, here we incorporate both (motor unit) temporal and spatial information to compute heterogeneous muscle activity. The computation of this motor unit information is described in Sect. 2.2.
The material-behaviour of the musculotendon complex also needs to account for connective tissues such as tendons and aponeuroses. This is achieved by generalising the muscle's constitutive relation with a tissue parameter where varies between 0 for pure tendon and 1 for pure muscle. To capture the different material responses of the tissues, the material parameters are linearly interpolated between values characterising pure tendon c t i and pure mus-  (4)

Enriching muscle activity with motor unit information
By assuming that muscle activity m (Eq. 4) can be fully characterised by the activity and distribution of muscle fibres at the microscale, we decompose the activation parameter accordingly (e.g. Röhrle et al. 2019), That is, we assume that the properties of all muscle fibres belonging to the same motor unit are identical and that they are activated synchronously. Then, individual motor unit activity can be described by a single, normalised, time-dependent function: where N MU is the total number of motor unit). We further assume that the spatial distribution of a motor unit in the muscle can be described by a volume fraction ̄i(x) ∈ [0, 1] ( x ∈ B ). Then, the activation of the muscle (at any point x and time t) can be computed by the ensemble of motor unit activity and distribution: Both motor unit activities and distributions are biophysically and micromechanically derived and are discussed in Sects. 2.2.1 and 2.2.2, respectively. By decomposing muscle activity in this way, we assume that the entire region occupied by the motor unit is activated simultaneously, i.e. the influence of propagating action potential is neglected. We investigate the error in motor output this introduces by comparison to a multi-scale, multi-physics model in Sect. 2.3.2.

Motor unit anatomy via virtual innervation
The virtual innervation process can be summarised as follows: (Step 1) reconstruct fibre microstructure from the fibre orientation field. (Step 2) define the number of motor unit and their size. (Step 3) perform virtual innervation iteratively to group the discrete fibres into motor unit. (Step 4) homogenise the innervated fibres within the finite element mesh.
(Step 1) Muscle architecture is represented via a fibre orientation field a 0 (X) ( X ∈ B ). We reconstruct idealised representations of the muscle microstructure ("fibre scaffolds") using an explicit streamline tracking method (e.g. Kupczik et al. 2015). This ensures that the reconstructed fibre scaffolds F k ( k = 1, … , N SC , where N SC is the total number of scaffolds) are compatible with the idealised continuous description of the muscle fibre architecture, i.e. the vector field describing the muscle fibre direction. A discrete representation of the scaffolds, i.e. given by a set of points X k micro ∈ F k , is iteratively computed by where r parametrises the scaffold, a 0 (X k micro (r)) is the fibreorientation at the previous scaffold coordinate, and Δ f is the step size along the fibre scaffold. Fibre scaffold tracking is initiated from a seed point s k ( k = 1, … , N SC ). The seed points are ensured to lie along the middle of the muscle geometry (retroactively). The tracking process yields a set of scaffolds F k with a corresponding set of seed points: S = {s k }. ( Step 2) The number of motor units to assign are specified by N MU . Motor unit sizes are characterised by their innervation ratio, which is the number of fibre scaffolds they contain. An exponential distribution of innervation ratios is used, based on Enoka and Fuglevand (2001), i.e.
Here, r IR = IR N MU ∕IR 1 is the ratio between the innervation ratio of the largest and smallest motor unit, IR tot = ∑ N MU i IR i is the total number of fibres in the muscle, and ⌈⋅⌉ is the ceiling function. ( Step 3) The virtual innervation can be considered as gathering N SC fibre scaffolds into i = {1, … , N MU } groups, with the number of scaffolds per group defined by IR s,i . Here, we introduce the subscript (⋅) i , which denotes that a quantity belongs to motor unit i, e.g. F k → F k i . This grouping process, or virtual innervation, proceeds iteratively from the smallest to the largest motor unit, i.e. first IR s,1 fibres are grouped, then IR s,2 , and so on, until IR s,N MU . A fibre scaffold (or seed point) can only be innervated by a single motor unit. First, the axon of a given motor unit i innervates a random seed point and this is set as the so-called reference seed point ( s k → s k i = c i ). Second, the remaining ( IR s,i − 1 ) axons continue to (randomly) innervate "free" seed points about the reference seed point c i , within a distance R, i.e.
where f d is a distance metric (discussed below) and S u is introduced as the set of currently uninnervated seed points. Once the desired number of seed points IR s,i have been innervated, the axon for motor neuron i + 1 innervates another (uninnervated) seed point and a new reference seed point is set ( c i+1 ). Then, the remaining axons of motor neuron i + 1 subsequently seek out free seed points about c i+1 (according to Eq. 10). The process continues for all remaining motor neurons. When a seed point is innervated, i.e. s k → s k i , the corresponding fibre scaffold including its microstructural points are also innervated, i.e. X k micro → X k micro,i . The search space for the reference seed point selection is restricted to a subregion of the muscle's cross section. This restriction is achieved by setting a seed point as an anchor point c a , and limiting the reference seed point search within a distance D about the anchor point, i.e.
This is analogous to Eq. 10. However, while R restricts the innervation of the axons of a given motor unit (about the reference seed point), D restricts the innervation of the reference seed points themselves (about the anchor point). The anchor point c a remains fixed for the entire innervation process. Note that it may be the case that as the innervation proceeds, no free seed points can be found within distances R and D. When this occurs, the selection distances R and D are incremented by ΔR and ΔD , and the searches are repeated.
The distance measure used to perform the search for (free) seed points is the generalised ellipsoidal distance, i.e.
where W is a weighting matrix which acts to scale the distance measure along its eigenvectors w 1 , w 2 , w 3 by the amount given by its eigenvalues, i.e. by 1 , 2 , 3 . For example, W = I (identity matrix) gives the Euclidean distance. In summary, the parameters {c a , R, D, ΔR, ΔD, W} characterise the virtual innervation procedure. ( Step 4) Thus far, motor unit anatomy has been described via discrete microstructural points. To map the innervated microstructure to the continuum-mechanical model, a statistical approach is used (e.g. Beran 1968;Kröner 1972). For a point in the muscle P ∈ B 0 with the coordinates X , a nearest neighbour search is used to group all fibre scaffold microstructural points X k micro,i (k = 1, … , N SC ) into (virtual) tissue samples. These samples are considered as statistical volume elements. Then, motor unit volume fractions ̄i(X) (i = 1, … , N MU ) are computed as the probability of finding a fibre scaffold point associated with motor unit i in the given statistical volume element, i.e.
where N micro is the number of microstructural points that have X as their nearest neighbour, of those N micro,i are the number of microstructural points that belong to motor unit i. Note that the motor unit volumes are "etched" into the muscle and are transformed to the current configuration via the placement function, i.e. ̄i(X) =̄i( (X, t)).

Motor unit activity
With the motor unit distributions at hand, we compute their activity for the duration of the simulation, t ∈ [0, T] . We assume that contractile properties and activity of a motor unit's fibres are identical. Therefore, a representative muscle fibre is used to characterise each motor unit's activity level. Recruitment and excitation-contraction coupling of the representative fibres is summarised as follows: (Step 1) A phenomenological model (Fuglevand et al. 1993) predicts motor neuron discharge times in response to volitional muscle excitation E(t). (Step 2) These discharges are coupled with a phenomenological membrane model (Aliev and Panfilov 1996) and (Step 3) a cross-bridge dynamics model (Razumova et al. 1999) simulates excitation-contraction coupling in the representative muscle fibre (Heidlauf et al. 2017, see also). This yields motor unit activities i (t) , which are unidirectionally coupled to the continuum-mechanical skeletal muscle model. That is, i (t) enters Eq. 7 and together with the motor unit fraction ̄i(X) determines the local muscle activation m (t, x) , which, in turn, enters Eq. 4. In the following, we only briefly summarise the recruitment-, calcium-and cross-bridge models; further details can be found in the literature, i.e. (Fuglevand et al. 1993;Razumova et al. 1999;Aliev and Panfilov 1996;Heidlauf et al. 2017). ( Step 1) The motor neuron discharge times t i,p , i.e. denoting the pth firing time of motor neuron i, are computed by supplying a global excitatory drive E(t) to the motor neuronpool. A motor neuron fires only when E(t) has exceeded its recruitment-threshold R i . The recruitment-thresholds are taken as proportional to the size of the motor unit (in terms of its innervation ratio) according to the size principle (Henneman 1957) (Fuglevand et al. 1993, c.f.). The mean frequency at which a motor neuron i fires is determined by where f p i and f m are the peak and minimum firing rates in Hz, respectively. ( Step 2) The discharge of a motor neuron gives rise to an action potential in the muscle fibre, which triggers the release of calcium (Ca 2+ ) into the sarcoplasm. This leads to cellular force production via cross-bridge cycling. In the phenomenological two-state membrane model of Aliev and Panfilov (1996), the fast-state variable represents the transmembrane potential. Further, it is assumed that the model's slow variable behaves qualitatively similar to the intracellular Ca 2+ concentration (e.g. Heidlauf et al. 2017). ( Step 3) Controlled by the Ca 2+ concentration, a simplified Huxley-type model (Razumova et al. 1999) simulates the muscle fibre's force output. Briefly, the model represents cross-bridges in four distinct states; with state transitions modelled by first-order kinetics. Although not described here, the rate coefficients dictate the twitch properties of the motor unit (or representative fibre). Here, rather than distinct categories of slow-and fast-twitch motor unit, we linearly interpolated the rate coefficients to create a smooth transition between the slowest-and fastest-twitch motor unit. Lastly, motor unit (or representative fibre) activity is assumed to be represented by the fraction of cross-bridge in the post-power stroke state ( A 2,i ), i.e. (14) where A max 2,i is the maximum amount of post-power stroke cross-bridge in the representative fibre (of motor unit) i.

Simulation trials
We used idealised geometries to investigate modelling assumptions prior to applying the method to an anatomically realistic musculoskeletal model as a proof-of-concept. Muscle activity is mesh dependent in the finite element mesh since it is computed by statistical averaging of a discrete representation of motor unit anatomy. We investigated this mesh dependency by using an idealised, 2D microstructural muscle model to compare "ground truth" activity (for various motor unit distributions) to that averaged at coarser mesh sizes (Sect. 2.3.1). Next, we analysed the impact of assuming instantaneous action potential propagation on force production of an idealised, 3D musculotendon complex by comparison to a multi-scale and multi-physics chemo-electro-mechanical skeletal muscle model (Sect. 2.3.2). Lastly, we used an anatomically realistic model of the masticatory system as a proof-of-concept to investigate the relationship between motor unit anatomy and bite force (Sect. 2.3.3)

Idealised geometry: mesh dependency of muscle activity
Motor unit volume fractions ̄i(X) are computed via statistical averaging within the finite element mesh (Eq. 13).
As a result, motor units are no longer spatially resolved at the individual fibre level but rather are represented by continuous distributions that are approximated with the finite element method. Since muscle activity is determined by motor unit distribution and activity (i.e. Eq. 7), the solution depends on the discretisation of the finite element mesh. We investigated this mesh dependency via an idealised, 2D, muscle cross section (Fig. 1). The analysis is summarised as follows: (Step 1) generate an idealised model of muscle cross section with muscle fibres (micrometre scale). ( Step 2) Virtually innervate the fibres to form motor units with different distributions. ( Step 3) Define various element sizes and compute the muscle activity within each mesh element. The solution for the finest macroscopic discretisation, termed as the reference element (RE), is used as reference to compute the error of the muscle activity.
(Step 1) An idealised virtual muscle cross section L CSA x = L CSA y = L CSA = 40 mm (cross-sectional area 16 cm 2 , Fig. 1 left) was populated with evenly spaced fibres at a distance of L fib x = L fib y = L fib = 100 m (Feher 2012;Polgar et al. 1973), giving a total of 16,000 fibres. Additionally, a reference element (RE) denoting the finest discretisation was defined as containing 4-by-4 fibres, yielding a RE size of Sharafi and Blemker 2011;Virgilio et al. 2015;Spyrou et al. 2017Spyrou et al. , 2019, c.f.) and totalling N RE = 10000 REs (Fig. 1 middle). ( Step 2) Innervation of the virtual fibres was carried out according to the algorithm described in Sect. 2.2.1. The motor unit pool was characterised by N MU = 50 , IR 1 = 100 , and IR tot = 10000 . Three different motor unit distributions were considered: completely random (maximal overlap), medium overlap, and minimal overlap. The respective shape parameters were: R = D = ΔR = ΔD equalling 2 L CSA , 0.25 L CSA , and 0.1 L CSA . The parameters c a = [L CSA ∕2, L CSA ∕2] and W = I were kept identical for all 3 distributions. Fig. 1 Ideal muscle geometry used for mesh dependency analysis. Each cross represents a virtual fibre, and the inter-fibre distance (or fibre diameter) is 100 µm. The reference element (RE) contains 4-by-4 fibres. Mesh sizes greater than RE sizes are considered, with an exemplary mesh size shown (Step 3) Statistical averaging was used to compute the motor unit volume fractions within each element for the different mesh sizes. With the volume fractions at hand, the muscle activity can be computed via Eq. 7. Different mesh sizes were defined in terms of the number of elements along a given side. Then the error in muscle activity was taken as the mean difference between RE and mesh activities, i.e.
where k = 1, … , N RE and the element e is found by a nearest neighbour search from RE k. Note that this error computation was repeated for each of the different mesh sizes considered. The mean per-mesh error was compared to mesh size for each of the 3 motor unit distributions considered. 1

Idealised geometry: modelling assumptions
The details of the multi-scale and multi-physics model (CEM model) used to investigate the influence of action conduction velocity can be found in Röhrle et al. (2012); Heidlauf and Röhrle (2014). Briefly, the CEM model consists of 1D muscle fibres embedded in a continuum-mechanical skeletal muscle model, where the fibres act to transmit the action potential from the neuromuscular junction outwards. We used the model to compare two cases, one where action potentials propagate along the fibres, and another, where the action potentials arise instantly along the fibres. The latter case reflects the modelling assumption used in the current model, i.e. synchronised activity along and between all fibres of a motor unit. The simulations were performed on a single core (Intel Core i7-4790K @ 4 GHz) with 24 GB of memory.
An idealised cuboid geometry with dimensions L × W × H of 40cm × 10cm × 10cm of a musculotendon complex was used. Both muscle and tendon tissue had the same cross section, with a muscle-to-tendon ratio of 2.5. The muscle was populated with 16 fibres, which were grouped in 5 motor units. The innervation ratios of motor unit 1-5 were 2;2;2;4;5, respectively. Virtual innervation was carried out with parameters D = R = ΔR = ΔD = W∕4 , with W = I and c a = [L∕2, W∕2, H∕2].
Stimulation of the motor units followed fixed firing rates of 28;17;14;10;7Hz for motor units 1-5, respectively. The motor units were recruited sequentially with a 10 ms offset from the preceding motor unit. For the first case, a current pulse was injected into the central node of the respective fibres. In the second case, an identical stimulus was delivered to all nodes of the respective fibres simultaneously. We adapted all parameters, besides the geometry, from the baseline experiment in Schmid et al. (2019). To compare the mechanical response of the two cases, the time taken for the force to attain a given value, defined as the time to force, was computed. Furthermore, electromechanical delay is calculated for both conditions as the time difference between the (first) stimulus and the time at which the musculotendon complex force increases by 1% of the maximum (twitch) force (Mörl et al. 2012).

Anatomical geometry: bite force specificity of masseters
The masticatory system model consists of the mandible, associated dental structures, and the left and right masseters. We systematically varied the motor unit territory distributions to produce multiple prototype models-each with unique masseter motor unit distributions. Then, each model was supplied with the same motor unit activity, and the resulting molar bite forces were taken as markers of motor output. Additionally, we simulated a static bite following the status-quo approach by applying an (equivalent) spatially constant activity to the masseters. The setup of the simulation is as follows: (Step 1) define masticatory model geometry and masseter motor unit properties. (Step 2) define the motor unit recruitment protocol. (Step 3) perform the static bite simulation. Lastly, motor unit anatomy was quantified via metrics outlined in Sect. 2.3.4. ( Step 1) Details of the masticatory model geometry and material properties are given in our previous publications (Röhrle and Pullan 2007;Saini et al. 2020;Saini and Röhrle 2022). Here, we focus on the definition of the masseter's motor unit territory. A total of approximately N SC = 1200 fibre scaffolds were reconstructed within each masseter. These fibre scaffolds were grouped into N MU = 50 motor units, with the smallest motor unit containing IR s,1 = 8 scaffolds and the largest IR s,50 = 53 (Eq. 9). That is, the formation of motor units 1 and 50 requires the innervation of 8 and 53 fibre scaffolds in the masseter, respectively.
All distributions shared two common features: First, smaller motor units were preferentially located anteriorly in the deep-head of the masseter by the placement of c a (Eq. 11) (e.g. Tonndorf and Hannam 1994;Eriksson and Thornell 1983). Second, the shape parameters 1 = 1, 2 = 3 = 5 were chosen such that ellipsoidal-shaped territories were formed with their main-axis aligned along the masseters' posterior-anterior axis (Eq. 12) (e.g. Tonndorf and Hannam 1994;McMillan and Hannam 1991) The remaining distribution parameters {R, D, ΔR, ΔD} were varied to produce n = 7 masticatory models (Table 1).
Lastly, the twitch properties of the motor units were as follows: the time to peak and peak twitch force for the smallest motor unit (number 1) were 0.091 N and 123 ms, respectively. The largest motor unit (number 50) had a time to peak a factor 3 smaller and a peak twitch force a factor 35 larger than the smallest motor unit. The properties for the remaining motor units were interpolated between the two extremes. ( Step 2) The masseters were recruited via a ramped excitatory drive E(t). It was linearly increased from 0 to 0.375 over Δt = 0.25s , held constant for Δt = 0.1 s , and then decreased linearly to 0 over Δt = 0.25 s . At its peak, the excitatory drive recruited 34 out of 50 of the motor units in the pool (selected motor unit activities are shown in 2a). The resulting motor neuron firing rates were characterised by Eq. 14. The maximum motor unit firing rate and range of firing rates ( f p 1 and f D , respectively) were 34Hz to 64 Hz (Kwa et al. 1995(Kwa et al. , 2002, the minimum firing rate f m was taken as 7Hz (Effect of motor 1989). The recruitment thresholds were adjusted such that approximately 50% of motor units were recruited for a 20% force level (Scutter and Türker 1998). Lastly, we computed an equivalent, spatially constant activity by a weighted sum of individual motor unit activity and the total normalised motor unit volume fraction of distribution n = 1 (Fig. 2b). ( Step 3) To simulate a static bite, i.e. contraction of the masseters does not displace the mandible, fixation boundary conditions were applied to the molar occlusal surface and the left and right fossas. The articular discs were fixated relative to the mandibular condyles and assumed to be in frictionless sliding contact with the fossas. The masseters were attached to the mandible over manually defined attachment areas and were fixated at their origins. The attachment areas extended caudally from approximately the middle of the ramus to the angle of the mandible. The molar was fixated over its occlusal surface, over which bite force was measured.
We recorded the bite force magnitude along the left-right, posterior-anterior and caudal-cranial axes at the right molar. Additionally, the angles made by the bite force vector with the horizontal and vertical planes were recorded. Simulations were performed using 32 cores (AMD Opteron-6373 @ 2.3 GHz) with 24 GB memory.

Motor unit distribution metrics
To quantify each n = 1, … , 7 motor unit territory distribution, three metrics were computed: the Motor Unit Neighbour Index (MUNI), the Motor Unit Territory Area (MUTA), and the mean active territory distance from the muscle centre of mass ( d c ). Note that these metrics were computed prior to homogenisation within the finite element mesh and representative metrics were taken from the right masseter. We computed the MUNI as follows: five random virtual tissue samples (dimensions: 5 mm × 5 mm ) were taken from the right masseter cross section, and the number of unique motor units within each tissue sample were  counted. These values were then averaged to compute MUNI for each distribution. We computed MUTA for each motor unit, i.e. MUTA n i ( n = 1, … , 7 , i = 1, … , N MU ), as the area of the convex hull over the fibres of a motor unit territory divided by the total cross-sectional area of the masseter (626 mm 2 ). Lastly, d n c was computed as the distance between the barycentre of the active motor units and the masseter centre of mass at the initial time step.

Mesh dependency of muscle activity
The mean error in muscle activity (Eq. 16) for the 3 motor unit distributions (maximum, medium, and minimal overlap) is plotted against mesh size in Fig. 3. At the coarsest mesh size, i.e. 1 element, the mean errors were 0.062, 0.13, and 0.22 as motor unit overlap decreased. Conversely, for the finest mesh size, i.e. when the element size is close to the reference element (RE) size, the means errors were 0.055, 0.047, and 0.030 as motor unit overlap decreased. The degree of mesh dependency was inversely proportional to the degree of overlap.

Influence of action potential propagation on force production
Computational times for the case with and without action potential propagation were 12 h and 6 h, respectively. The reaction forces for both cases, summed over the nodes at the fixed tendon end, are shown in Fig. 4. The force response for the action potential propagation case is slightly smoother and delayed in time compared to the non-propagation case. The difference in time to force between both models was 16.6ms. Both models reached the same maximum force magnitude of 39N. The electromechanical delay for the cases with and without action potential were 18.8 ms and 8.2 ms, respectively.

Masseter motor unit anatomy
Examples of the virtual neuromuscular anatomy are shown over the masseter cross section in Fig. 5. Further, the quantitative motor unit metrics for all motor unit distributions are given in Table 2. For small reference seed point search parameters {D, ΔD} , motor unit types were confined to subregions of the masseter cross section. For example, distributions 1, 4 and 5 ( {D, ΔD} = {10 mm, 2 mm} ) show type-S (smaller) motor units positioned anterio-medially in the masseter (Fig. 5). On the other hand, larger values {D, ΔD} ={300 mm,2 mm} resulted in scattered motor unit territory (distribution 6 in the aforementioned figure).
For larger motor unit axon search parameters {R, ΔR} , more overlap was observed between neighbouring motor unit territory. For example, distribution 4 shows a larger overlap than distribution 1 (Fig. 5). This is also reflected in the motor unit territory features.  Fig. 3 Mesh dependency of the spatial distribution of muscle activity in the idealised mesh. The error is the mean difference between muscle activity at the reference element (RE) and mesh elements. The three curves represent different degrees of motor unit overlap; circular marks for maximum overlap (random distribution), plus marks for medium overlap, and triangle marks for minimal overlap. The vertical dotted lines show the possible extremes of mesh size Fig. 4 Tendon force of the idealised geometry for the model with action potential propagation (dotted, red) and without action potential propagation (solid, black). The motor unit firing times are indicated by the vertical ticks for each motor unit (right axis) corresponded to {R, ΔR} = {1 mm, 4 mm} , yielding more segregated motor units. The reference seed point search parameters {D, ΔD} influenced MUNI to a lesser degree. A similar relationship between the distribution parameters and MUTA was observed. That is, the axon search parameters {R, ΔR} correlated (positively) to a much higher degree to MUTA than did {D, ΔD}.
Note that for the spatially constant activity case, each of the 50 motor units exists at all points. Therefore, MUNI const = 50 , and since each motor unit territory covers the entire muscle, MUTA const for each territory and the cross-sectional muscle area are identical; at 100 %.

Bite force specificity of masseters
Computational times for all static bite simulations were between 10 h and 15 h. Cross-sectional masseter activity at point of maximum force production is shown for selected distributions in Fig. 7. The cross sections reveal heterogeneous activity, which differs among the distributions e.g. distributions 1, 4 and 5. The differences in masseter activity patterns also produce different bite forces. The horizontal bite force component at the corresponding time instance is overlaid in Fig. 7.
The distance between the barycentre of contraction and the masseter centre of mass d c is given in Table 2. Smaller {D, ΔD} resulted in larger d c values (4.4-6.1 mm) (distributions 1,2 and 4), meaning the barycentre of contraction was positioned further away from the masseter centre of mass. Conversely, as {D, ΔD} increased, d c tended towards zero, e.g. distributions 3, 5 and 6 all have d c < 1 mm.
The evolution of the bite force components along the left-right, posterior-anterior and caudal-cranial axes, respectively: F n LR , F n PA and F n CC , for n = 1, … , 7 distributions, are plotted in Fig. 6, together with the force components of the spatially constant activity model. The peak force magnitude for the spatially constant activity model was F const M = 147N . The peak force magnitudes for the n = 1, … , 7 distributions had a mean and range of 146 N and 5 N, respectively. The peak forces in the left-right, posterior-anterior and caudalcranial axes had means of 14N, 1 N, and 145 N, respectively. The corresponding ranges were 6N, 4 N, and 6 N.
The evolution of the angle of the horizontal bite force component for the n = 1, … , 7 distributions is shown in Fig. 6b, together with the horizontal angle for the spatially constant activity model. At the point of maximum force, the horizontal angle for the spatially constant activity model was const M = −70 • (note that 0 • corresponds to the anterior direction and − 180 • to the posterior direction). The horizontal angles for the n = 1, … , 7 distributions had a mean and range of − 86 • -40 • , respectively. For the spatially constant activity model, the force was oriented anteriorly and showed minimal variation during the (ramp) contraction. In general, the motor unit driven models produced bite-forces which were less steady in their horizontal angle, with some models even showing a reversal in bite-force direction along the posterior-anterior axis, i.e. distributions 5 and 6. Bite forces of distributions 1 and 2 were oriented posteriorly during the entire contraction, and distributions 3, 7, and 4 were oriented anteriorly during the entire contraction.

Discussions
In some regards, continuum-mechanical models of skeletal muscles (3D models) assume a lumped function. For example, a (single) representative muscle fibre is used to characterise the entire muscle's activity level. However, muscle contractile properties and activity are heterogeneous due to the muscle's organisation in motor units. This means that, the selective recruitment of regionally confined motor units can allow task specific muscle activations. Here, we presented a method to introduce such heterogeneity within 3D models by incorporating motor unit distribution and activity information. While this incurred minimal additional computational costs, it did increase Table 2 Force magnitude ( F M ) and horizontal orientation ( H ) at point of peak force production. Motor unit territory metrics: Motor Unit Neighbour Index (MUNI), Motor Unit Territory Area (MUTA) (masseter cross-sectional area 626 mm 2 ), and centre of contraction d c are also shown for each distribution. The subscript 10 and 90 refer to the 10th and 90th percentile of the MUTAs. Distribution "const." refers to the status quo, spatially constant activation model model complexity. The question then arises: when is the inclusion of motor unit information in 3D models paramount to muscle motor output predictions? We address this question in regard to motor unit distribution and functional heterogeneity of muscles. Having shown the conditions under which such a motor unit enriched 3D model is beneficial, we conclude with the limitations of the current model and its implications for continuum-mechanical muscle models.

Motor unit distribution and bite force
Muscles contain hundreds to thousands of motor units, and motor units, in turn, contain tens to thousands of muscle fibres, and represent the smallest unit that can be controlled individually by the nervous system. A motor unit's fibres are often regionally confined within the muscle. For example, higher numbers of small motor units were found deep in the masseter (Eriksson and Thornell 1983;Tonndorf and Hannam 1994;van Eijden and Turkawski 2001), or superficially in the tibialis anterior (Mesin et al. 2010). The sequential recruitment of locally confined motor units results in localised muscle contraction. Accordingly, fibre orientation within this region dictates the force vector produced. For example, Turkawski et al. (1998) found, by individually recruiting motor units, that the rabbit masseter produced "lines of action ... at least as large as the variation in fibre directions". Other studies show a similar link between regional muscle (or motor unit) activity and changes in force vector (Schindler et al. 2005;Ogawa et al. 2006;Holtermann et al. 2009;Miyamoto et al. 2012). This work presents a computational model to demonstrate the relationship between motor unit anatomy, recruitment, and bite force using anatomically realistic masticatory model. The results showed that bite force magnitude was mildly affected upon variation of recruitment regions (range of 5 N, for an average 146 N); however, significant differences were observed in bite force direction (range of 40 • , for an average of − 86 • ). Furthermore, differences in force direction seemed to increase as the motor unit territories were less overlapped. For example, distributions 1 and 5 had highly segregated motor units (motor unit neighbour index (MUNI) of 6.7 and 7.1, respectively) and produced a highly oscillating force direction compared to distributions 4 and 2 (MUNI of 18.4 and 16.6, respectively) ( Table 2 and Fig. 5).
Recalling that high MUNI values correspond to increased motor unit overlap. As motor units become more overlapped, their recruitment results in the excitation of the same muscle region, therefore producing a more stable force direction. Indeed, as MUNI values increased, from 6.7 to 15.6 to 18.4 (distributions 1, 2 and 4, respectively), the steadiness of bite force direction increased (Fig. 6b). Taken to the extreme, all motor units would overlap at every point within the muscle; this is the assumption behind the spatially constant activity case (which has a theoretical MUNI of 50). As expected, this produces the most steady force direction (Fig. 6b).
The link between motor unit overlap and force direction reveals the interplay of muscle architecture and neural organisation. That is, a muscle's force direction can only be modulated by the nervous system when fibres of the motor units are regionally confined. The masseter has more confined territories as compared to muscles of the limbs -"suggest[ing] a more localised organisation of motor control in masticatory muscles" ( van Eijden and Turkawski 2001). Additionally, Schindler et al. (2014) found that small "sub-volumes" of the masseter were recruited under particular tasks via the recruitment of adjacent motor units. This suggests that modelling masseter activity (and that of other muscles with complex architecture) as spatially constant may oversimplify force predictions.
In these static bite simulations, masseter architecture explains the relationship between motor unit overlap/ segregation and force (direction) steadiness. The masseters were defined as having a smaller deep head and a much larger superficial head. The deep (or inferior) head arises from the deep surface of the zygomatic arch and inserts on the upper part of the ramus. The superficial (or superior) head is much thicker and arises from the anterior two-thirds of the zygomatic arch and inserts to the angle of the mandible. When viewed from the sagittal plane, fibres of the deep and superficial heads form an "x", i.e. they are aligned posterio-cranially and anterio-cranially, respectively. In the bite simulation, the same set of (34 out of 50) motor units were recruited via identical firing patterns. The only difference was the distribution of their fibres within the masseters (Fig. 5, bottom row). When coupled with the different fibre directions in the deep and superficial head, this means that the relative location of the contracting fibres influences the line-of-action. Figure 7 shows the pattern of activity over cross sections of the left and right masseters, at the point of maximum recruitment -for different motor unit distributions. The activity patterns correspond well with the motor unit fibre distributions in Fig. 5 (bottom row), as expected. The horizontal bite force angle is also shown in Fig. 7, and reveals that as the recruited fibres become more dispersed throughout the masseter (from distribution 1 to 6), the force direction converges to the spatially constant activation case. More specifically, distributions 1 and 2 had contracting motor units in the deep head, and thus the bite force was directed posteriorly (horizontal angles of − 114 • and − 100 • , Table 2). Note that a horizontal angle of 0 • , − 90 • , and − 180 • refers to forces directed anteriorly, to the left, and posteriorly, respectively (Fig. 7). Distributions 3, 5, 6, and 7 produced anteriorly directed bite forces (horizontal angles between 70 • and 80 • ) since most of the contracting motor units were located in the superficial head. Interestingly, distributions 5 and 6 produced a bite force that switched direction; posteriorly to anteriorly. Here, motor units are initially recruited from the deep head and later from those located in the superficial head. These findings are in agreement with experimental observations, which show the involvement of the masseter's deep and superficial regions to retract and protrude the mandible, respectively. For example, Belser and Hannam (1986) found that during jaw retraction, the respective activity of superior and deep fibres was 5.5 % and 47.5 %. Additionally, Hannam and McMillan (1994); Blanksma et al. (1997) found that deep fibres in the masseter contributed predominantly to jaw elevation and jaw retraction. Recording the activity of individual motor units, Ogawa et al. (2006) found that those motor units located deep within the masseter started firing earlier when the bite was oriented "posterio-laterally", i.e. retraction of the jaw laterally. As the bite force became more anterior, the motor units in the superficial region became more active.
Although distributions 5 and 6 produce similarly oriented bite forces, the motor unit territory in distribution 6 show much greater overlap than distribution 5. This implies that motor unit territory overlap alone does not characterise the force behaviour; rather the position of the motor unit territory is of paramount importance. This is reflected by the fact that distributions 5 and 6 have a similar mean barycentre of contraction ( d c in Table 2), again, despite having differing amounts of territory overlap. Conversely, distributions 2 and 7 have similar amounts of overlap and territory areas; however, since they differ in their mean barycentre of contraction, one produces a slightly higher force, oriented posterio-cranially (distribution 2), while the other is oriented anterio-cranially. Subtleties of muscle function as discussed above cannot be captured without taking motor unit distributions into account.

Assumptions and limitations
Force production within a fibre follows the propagation of the action potential. As with status-quo 3D models, we assume that action potential propagation time (in the order of milliseconds) is negligible when compared to the temporal scale of movements (in the order of seconds). When action potential propagation was ignored, the force detected at the tendons was advanced by approximately 20ms. The geometry and material properties of the musculotendon complex impact this value. For example, Schmid et al. (2019) showed that a low muscle-to-tendon ratio, long fibres, stiff muscles and soft tendons act to (individually) increase electromechanical delay. Therefore, care should be taken when action potential propagation is neglected for muscles fitting this profile.
Muscles produce force by the orderly and sequential recruitment of their motor units. Motor unit recruitment may, however, be altered in some cases, such as the task performed, contraction velocity, and sensory feedback (see Hodson-Tole and Wakeling 2009, for example). Alternate recruitment strategies may enhance regional contractions since the nervous system could recruit neighbouring motor units. For example, Schindler et al. (2014) found localised task groups of motor units in the masseter. In the current model, motor units are recruited strictly according to the size principle, i.e. modification of recruitment thresholds is not accounted for.
The homogenised continuum-mechanical model does not explicitly resolve the muscle fibres. Instead, the most important properties of a muscle's fibre architecture, in terms of mechanical force production, are included in our model in a continuous manner. Specifically, the fibre orientation field and the motor unit volume fractions. Similar assumptions have been used in the past for chemo-electro-mechanical, multi-physics models of skeletal muscle tissue. This is mostly due to the fact that models explicitly resolving the full microscopic muscle geometry are typically limited to simulate small tissue samples due to computational costs (Sharafi and Blemker 2011;Virgilio et al. 2015;Spyrou et al. 2017Spyrou et al. , 2019. However, a continuous representation of motor unit anatomy introduces mesh dependent model error. The mesh dependency study revealed that the mean error in muscle activity depends on both mesh size and motor unit overlap. For maximal amount of overlap, i.e. randomly distributed motor units, the error did not show convergence with mesh refinement. This is because each element has a very high probability of containing all the muscle's motor units. Therefore, as element size changes it has little to no effect on the mean activity within the element. On the other hand, minimally overlapped motor units show a strong reduction in error upon refinement. This is because low motor unit overlap leads to regions in the muscle that contain only a single motor unit. Then, as the mesh size decreases, the probability that an element covers a region that contains only a single motor unit increases. This explains the observation that the error for the minimally overlapped case drops below that of the random distribution case, since random motor unit distribution will not yield elements that contain a single motor unit. Regarding the proof-of-concept masticatory simulation, the element edge lengths in the masseter ranged between ≈ 2 − 3 mm , which corresponds to roughly 6 % mean error in the computed muscle activity. It should be noted that mesh size does not influence muscle activity alone, but rather also impacts the fibre orientation and stress-strain fields. A comprehensive convergence analysis is an important next step in verifying the proposed method, and finite element skeletal muscle models in general. We assumed that fibre size is the same within and between motor units. It has been observed, however, that fibres between motor unit types differ in size (e.g. Polgar et al. 1973) and have normally distributed diameters within a motor unit type (e.g. Hegarty and Hooper 1971). Furthermore, contractile properties also seem to vary between motor unit types. Kanda and Hashizume (1992) showed, first, that mean muscle fibre cross-sectional area increased proportionally with innervation ratio. Second, as innervation ratio increased, so did the specific tension (motor unit tension divided by motor unit cross-sectional area area). Therefore, not only do larger motor units occupy more space, but they produce a greater force per unit of area. Factors which govern the peak tension of a motor unit in the current model are activation, specific strength, and innervation ratio. The activity of each motor unit is normalised (Eq. 7) between [0, 1] and the specific strength is not specified per fibre, but rather by the peak isometric stress, i.e. S max , and is constant over the entire muscle (Eq. 4). Therefore, neglecting geometrical effects and assuming full activation of all motor units, differences in peak force between motor units depend solely on innervation ratio. While the inclusion of these motor unit characteristics would yield more accurate motor unit force production, it appears that the innervation ratio of a motor unit is the dominant factor in predicting its force output (Kanda and Hashizume 1992). Therefore, from a functional point of view, although fibre diameters and specific tension are idealised, the relative motor unit force contribution is maintained in the current model.

Implications
The interplay of motor unit anatomy, selective recruitment by the nervous system, and the architecture of the muscle gives rise to regional contraction and functional heterogeneity of the muscle. From a theoretical point of view, this would be most apparent for muscles with "complex" architecture (bipennate or multipennate) and those with a large number of segregated motor units. We hypothesise that the more versatile a muscle needs to be, the more complex its architecture and the more segregated its motor units. This is intuitive since a complex architecture can only be exploited by the nervous system if the motor units are regionally confined.
Joint kinematics typically result from the coordination of multiple spanning muscles. One way to determine a muscle's contribution to a movement is via inverse kinematics and optimisation techniques. Typically, modelling studies assume each muscle to operate as a single contractile unit, whether using 1D Hill-type models or 3D models (Péan et al. 2019;Harandi et al. 2020). However, by considering the joint as spanned by groups of motor units instead, the spatial resolution of force production would be greatly increased. Consequently, the sub-maximal contraction of motor units would no longer simply scale the resultant force vector's magnitude but would also influence its direction and thus the joint torques produced (as shown by the bite force in the present study). This has further implications for muscle function and coordination, such as fatigue during sustained contractions. Here, the coordination of motor units across multiple muscles could be optimised in such a way as to maintain the posture, e.g. via recruitment strategies such as motor unit substitution or rotation (Bawa and Murnaghan 2009).

Conclusions
In this paper, we developed a method to generate and include motor unit activity and anatomical information in 3D continuum-mechanical muscle models (3D models). Although idealised, the proof-of-concept simulation study showed that enriching 3D models with motor unit information significantly impacted bite force. In summary, (i) activating different regions of the masseter had a minimal effect on bite force magnitude, but affected bite force direction significantly; (ii) activation of the deep and superficial regions of the masseter yielded a more posteriorly and anteriorly directed bite force, respectively; (iii) lower overlap between motor units yielded a higher variation in bite force direction; and (iv) motor unit overlap has a smaller impact on bite force (magnitude and direction) than motor unit position. These findings suggest that the current modelling approach would benefit muscles with complex internal architecture and segregated motor units. The proposed method was shown to reasonably approximate skeletal muscle physiology. Nevertheless, a rigorous validation is still required. In particular, the in vivo characterisation of motor unit territories presents a unique challenge. Despite this, the present modelling framework enables, at least in a qualitative sense, the investigation of motor unit properties within the context of muscle motor output and joint function. This is relevant not only for understanding the functional heterogeneity of healthy muscle but also for aged or pathological muscle, where motor unit remodelling takes place.
Funding Open Access funding enabled and organized by Projekt DEAL. This work was supported by the German Research Foundation (DFG) project GRK2198 (Soft Tissue Robotics), and by the German Federal Ministry of Education and Research (BMBF) project 01EC1907B (3DFoot).

Conflict of interest
The authors have no competing interests to declare that are relevant to the content of this article.
Ethics approval Not applicable.

Consent for publication Not applicable.
Code availability Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.