SUMMARY
This paper reviews current mathematical models of sphincters and compares them with a new spatial neuromuscular control model based on known physiological properties. Almost all the sphincter models reviewed were constructed as a component of a more extensive model designed to mirror the overall behaviour of a larger system such as the lower urinary tract. This implied less detailed modelling of the sphincter component. It is concluded that current sphincter models are not suitable for mimicking detailed interactions between a neural controller and a sphincter. We therefore outline a new integrated model of the biomechanics and neural control of a sphincter. The muscle is represented as a lumped-mass model, providing the possibility of applying two- or three-dimensional modelling strategies. The neural network is a multi-compartment model that provides neural control signals at the level of action potentials.The integrated model was used to simulate a uniformly activated sphincter and a partially deficient innervation of the sphincter, resulting in a non-uniformly activated sphincter muscle. During the simulation, the pressure in the sphincter lumen was prescribed to increase sinusoidally to a value of 60 kPa. In the uniformly activated situation, the sphincter muscle remains closed, whereas the partially denervated sphincter is stretched open, although the muscle is intact.
Introduction
The main mechanical functions of sphincter muscles are controlled constrictions and relaxations. Sphincters are used for temporary mechanical separation of compartments ranging from a valve-like function in one-directional transport of food in the gastrointestinal tract to a water-tight seal of the urinary bladder or the rectum. Different functions and topologies demand architectural adaptations of a sphincter and adaptations of its neural control.
Mathematical models of sphincters can be used to study various aspects of these adaptations. In most current models, the sphincter is implicitly included in a larger system, for instance in a model of the lower urinary tract (Bastiaanssen et al., 1996a; Bastiaanssen et al., 1996b). In these systems, the model of the sphincter is a simple structure, because the complete models were designed to simulate the collective behaviour of a larger system and were not intended to study the detailed internal dynamics of a sphincter. Most published models of sphincters include only one ‘muscle unit’. The muscle units range from a resistance with an on/off switch control to a resistance with muscle-like dynamics and most have a one-dimensional geometry. To study detailed aspects of a distributed parameter system, such as a sphincter or a neural network, more complex spatial models of several interacting units are needed. We found only one spatial model of a sphincter in the literature (Gielen, 1998).
After reviewing current models of sphincter, we outline a new integrated neuromuscular model of a sphincter and its neural control that allows detailed study of various aspects of this system, such as the spatial distribution of different muscle fibre types and stretch receptors, the effects of spatial neuromuscular defects, neural activation and dynamics. This integrated neuromuscular model consists of two parts. The muscle is represented as a lumped-mass model offering the possibility of applying two- or three-dimensional modelling strategies. The neural network is a multi-compartment model that examines signal flow in the neural controller at the level of action potentials.
Comparison of models in literature
Sphincter muscle models
Current sphincter models are compared in Table 1. All current models of sphincters are portions of models designed to represent the lower urinary tract, except for a model of the pupil of the eye (Usui and Hirata, 1995) and the finite-element model of Gielen (1998). In many models of the lower urinary tract, the sphincter is not modelled separately but is integrated into the urethra (Drolet and Kunov, 1975; Fröhlich et al., 1977; Hosein and Griffiths, 1990; Valentini et al., 1992; Van Duyl, 1993; Hübener and Van Mastrigt, 1994). In these models, closing and opening of the urethra are realised by multiplying the outflow by the value of a simple block function (with possible function values of 0 or 1) (Fröhlich et al., 1977; Van Duyl, 1993) or by applying a simple pressure threshold that, when exceeded, allows urine outflow (Drolet and Kunov, 1975; Hosein and Griffiths, 1990; Valentini et al., 1992; Hübener and Van Mastrigt, 1994).
In some models (Hosein and Griffiths, 1990; Van Duyl, 1993; Hübener and Van Mastrigt, 1994), the mechanical behaviour of the urethra is represented by a non-linear relationship between bladder pressure and flow rate, known as the urethral resistance relationship (Griffiths, 1980).
In three models (Usui and Hirata, 1995; Bastiaanssen et al., 1996b; Van Duin et al., 1998; Van Duin et al., 1999), a different approach was followed. A sphincter was modelled using visco-elastic elements and activity-dependent contractile elements (with included visco-elasticity). These models show muscle-like dynamics and will be discussed in the next section. Gielen (1998) proposed a spatial finite-element model for skeletal muscles that was also used to represent a sphincter.
Single-unit models with muscle dynamics
Hill’s three-element visco-elastic model is a commonly used representation of muscle. The mechanical behaviour of the muscle is described by a Hill-type contractile element (CE) based on Hill’s characteristic equation (Hill, 1938), which was later shown to be a phenomenological description of the behaviour of the actin and myosin proteins (Huxley and Niedergerke, 1954; Gordon et al., 1966). It includes elastic elements (EEs) (see Fig. 1A). If the CE is not activated, the EE represents the passive elastic properties of the muscle.
Knowing the strain rate of the elastic elements and active state, the mechanical stress in the muscle can be derived (see Fig. 1A). This stress is assumed to be uniform in the muscle. The normal stress on the luminal fluid generated by the sphincter muscle can be calculated when its state of activation, radius and width are known. The normal stress is assumed to be in force equilibrium with the stress on the muscle as a result of the intraurethral pressure (Van Duin et al., 1998; Van Duin et al., 1999), which implies that inertial effects are ignored. This type of model is therefore known as an equilibrium model. Bastiaanssen et al. (1996a,b) also used an equilibrium model, with one EE (see Fig. 1B). The force production of the CE in the sphincter depends on the length of the CE, the specific shortening velocity and the state of activation of the CE. The passive EE and the active cross-bridges contribute, in parallel, to the tensile stress in the sphincter. A homogeneous distribution and shortening velocity of the CE and the EE in the muscular layer of the sphincter is assumed. In addition, a uniform activation of the CE is assumed. As a result, the tensile stress in the sphincter can be described by the stress production of a single muscle fibre. The tensile stress of the single muscle fibre equals the sum of the active tensile stress due to the cross-bridges in the muscle fibres and the passive tensile stress of connective tissue and epithelium due to the EE.
Usui and Hirata (1995) included a viscous element (VE) in the muscle unit of their model pupil (see Fig. 1C). The sphincter and dilator models (the pupil includes a sphincter and a dilator muscle with radially directed muscle fibres) were both constructed using the same basic units (see Fig. 1C), but using different parameters, and were combined according to the structure of the pupillary muscle. The pupillary muscle system was simplified to a one-dimensional push/pull structure using Newton’s second law to calculate the variation in pupil radius from the forces exerted by the sphincter and dilator muscles on a characteristic mass.
Multi-unit model with muscle dynamics
Gielen (1998) proposed a continuum description for a contracting skeletal muscle in which the distributed moments model of Zahalak (1981) was implemented to derive the constitutive equations. In this approach, the cross-bridge forces are essentially described using the Huxley (1957) description. The force of a half-sarcomere is derived by summation of the contribution of the available cross-bridges. The continuum description was solved numerically using the finite-element method. Mass inertia was ignored in this model. This approach is promising, but a test of the model (Gielen, 1998) with a prescribed axi-symmetric load and an initial axi-symmetric sphincter geometry was only partially successful because significant deviations were found from the ideal axi-symmetric solution. The deviations were attributed to limitations in the applied mesh and the associated choice of the properties of the elements.
Neural control of the model sphincter
The neural control systems of current sphincter models are compared in Table 1: some models lack a neural controller for the sphincter (Drolet and Kunov, 1975; Fröhlich et al., 1977; Van Duyl, 1993; Gielen, 1998), and the process of opening and closing depends only on an on/off switch that, in some models, was coupled to a threshold pressure (Fröhlich et al., 1977; Van Duyl, 1993). Valentini et al. (1992) and Hübener and Van Mastrigt (1994) prescribed the neural excitation of the urethra as a function of time. This approach is known as an open-loop control system.
Closed-loop neural control
Hosein and Griffiths (1990) postulated a neural control system involving two mutually inhibitory control regions, leading to a bi-stable (storage-voiding) system. The reflexes were driven by a single afferent signal with contributions from bladder wall tension and urethral distension.
Usui and Hirata (1995) modelled the portion of the autonomic nervous system innervating the pupil: the parasympathetic and sympathetic input to the sphincter and the dilator muscle, respectively. The parasympathetic and sympathetic nervous activities were mutually inhibitory.
Bastiaanssen et al. (1996a) used a qualitative model of the neuronal circuitry involved in the control of the uropoietic system (Kinder et al., 1995) to build a detailed quantitative model of the different anatomical structures involved in the system.
Van Duin et al. (1998, 1999) presented a less complicated alternative description of the neural control of the lower urinary tract. The input parameters to the neural model were related to the length of the bladder and urethral wall, to the pressure in the bladder and to the stress in the bladder wall.
Current neuromuscular sphincter models
All previous models were designed to mirror the overall behaviour of a large system. This implied less-detailed modelling of components so that, for instance, the pressure in the sphincter is reflected by only one value. For some purposes, this coarse modelling is appropriate. However, the pressure in a sphincter muscle follows a gradient, rather than having a single value, and the state of activation varies over the sphincter (Podnar and Vodušek, 1999; Podnar et al., 1999; Podnar and Vodušek, 2000). In the literature, the term distributed parameter system is often used for this type of system. Distributed parameter systems can be modelled by a continuum approach (e.g. Gielen, 1998). Alternatively, spatially discrete (lumped) models can be used in which continuously distributed parameters (properties) are concentrated in the crossings of a spatial grid (spatial discretisation). For both approaches, more details can be revealed by using a finer mesh or grid at the expense of higher computer memory requirements and a longer computation time. By expressing the state of a sphincter in a single value, much information is ignored. To study the finely distributed properties of a typical distributed parameter system, less coarse models are needed. Current sphincter models are not suitable for mimicking the detailed interactions between neural control and sphincter.
In the next section, we propose a spatial model for the neuromuscular control of sphincters in which a multi-compartmental model of a neural network is connected to a lumped-mass model of a sphincter muscle.
The lumped-mass sphincter model
The sphincter is modelled as a two-dimensional system of n point masses. The masses are assumed to be interconnected by tensile elements (with zero cross-sectional area), thought to represent muscle fibre bundles and connective tissue. The topology of masses and tensile elements is shown in Fig. 2A. The masses are connected to either four (corner point) or two (mid point) tensile elements, with the exception of the corner point masses on the inner or outer border, which are connected to only three tensile elements. The tensile elements are assumed to have a constant curvature along their length, as determined by the locations of two adjacent corner points and one intermediate mid point. Fluid elements are enclosed by the tensile elements through four corner points and four mid points. The fluid is assumed to be slightly compressible, described by a bulk modulus. The pressure is assumed to be constant in any one element. Both the active and the passive elastic elements have non-linear properties. For the active component of muscle tissue, non-linear force/velocity (Hill-type) relationships have been assumed, whereas length-dependence on filament overlap (and thus force development) is included. The stress of the passive elements is non-linearly related to the strain. Tensile elements are assumed to be arranged in a circular fashion. Muscle fibres may extend over one or more tensile elements. Their tensile stresses are assumed to depend on passive tensile properties, passive viscous properties and activity-related dependencies of longitudinal strain and strain rate, described after Van Leeuwen and Kier (1997). The activity of the tensile elements is described by first-order dynamics with the local axonal firing rate as input (see Bastiaanssen, 1996). The radial tensile elements are described by passive elastic and viscous properties only. The motion of the masses is calculated by application of Newton’s second law. The acceleration of each mass is calculated from the forces in the connecting tensile elements and the forces due to the pressure in the adjacent fluid compartments. Briefly, the set of n coupled second-order differential equations is first rewritten to an equivalent system of 2n first-order differential equations and then numerically solved using the Adams–Bashford–Moulton method (Burden et al., 1981; Mathews, 1992).
Neural network model
Compartmental modelling
When constructing detailed neural models that explicitly include many of the potential complexities of a cell, the usual approach is to divide the neuron into a finite number of interconnected compartments. Each compartment is modelled with equations describing an equivalent electrical circuit (Rall, 1959). With the appropriate differential equations for each compartment, depending on its active and passive electrical properties, one can model the behaviour of each compartment and its interactions with neighbouring compartments.
Conceptually, electrically short segments are assumed to be isopotential and are lumped into a single membrane compartment with a specified resistance and capacitance, either passive or active. Compartments are connected to one another via a longitudinal resistivity according to the topology of the tree describing the detailed morphology of the original neuron. Hence, differences in physical properties (e.g. diameter, membrane properties, etc.) and differences in potential occur between compartments rather than within them. The length of a compartment must be small enough for its electrotonic length to be greater than 0.2λ (Segev et al., 1985) to consider it as isopotential. The space constant λ is the length of a cable with the same diameter (as the compartment) that has a membrane resistance equal to the cytoplasmic resistance. Often, this requires simulations with very large numbers of compartments. A GENESIS (GEneral NEural SImulation System) (Wilson et al., 1989; Bhalla and Bower, 1992; Bower and Beeman, 1994) model of a cerebellar Purkinje cell was recently published that uses 4550 compartments and 8021 channels (De Schutter and Bower, 1994a; De Schutter and Bower, 1994b).
The computation time for a simulation depends on several factors: the integration method, the precision of the method, the size of the model and the total number of conductances. One can simplify the model of a neuron to keep the number of compartments (and the computation time) reasonable. The classical way to simplify dendritic trees is to lump branches together into an equivalent cylinder, i.e. to use one large dendrite instead of several small ones. An extreme of this method is the ‘ball-and-stick’ model used in linear cable models: one spherical somatic compartment and a large cylindrical compartment equivalent to the axon, the dendrites, etc. This ‘ball-and-stick’ method was used to construct the neuron model in the following section. The simplification of the motor neurons was motivated by the fact that a more complicated model would make no difference in the current simulation. In this simulation experiment, the motor neurons are responsible only for the uniform and synchronous activation of the sphincter muscle.
Neural structure
The model neuron consists of a somatic segment, a dendritic segment and one axonal compartment. The dimensions of the soma compartment were based on reconstructions of external urethral sphincter motor neurons from Onuf’s nucleus (Sasaki, 1994). In the present model, the axon is represented by a sphere with a diameter of 28.40 μm. Axonal segment length and diameter are 100 μm and 5 μm, respectively. The length and diameter of the dendritic segment are 100 μm and 4 μm, respectively.
Passive and active parameters
The passive parameters of the neurons were set according to the values reported by Fleshman et al. (1988) for cat α-motor neurons. Specific membrane capacitance Cm was 1.0 μF cm–2, and specific cytoplasmic resistance Ri was 70 Ω cm. The neuron had values for the specific membrane resistance Rm of 11 kΩ cm2 for the dendrite and 0.225 kΩ cm2 for the other compartments (soma and axonal segment).
Active properties were associated with the initial segment and soma compartment only and consisted of five voltage-dependent channel types modelled using Hodgkin–Huxley equations (Hodgkin and Huxley, 1952). The conductances in the neuron segments are given in Table 2. The axonal segment contained a fast Na+ conductance and a fast delayed rectifier (Barrett and Crill, 1980). The soma compartment contained a fast Na+ conductance, a fast K+ conductance and a slow K+ conductance. The K+ conductances are analogous to those described by Barrett et al. (1980) for cat motor neurons. Jones and Bawa (1997) used the same composition for a comparison between cat and human motor neurons. They suggested that the biophysical properties governing rhythmic firing in human motor neurons are similar to those of the cat (see Table 2 and Table 3 for a detailed description of parameters).
A postsynaptic conductance was implemented in the dendritic segments, mimicking the time delay caused by transmitter release and offering the possibility of simulating input from other neurons. The model also offers the possibility of injecting current in an arbitrary segment or neuron. The neuron model offers a choice between two different interfaces between the neural network and the muscle, a voltage-dependent presynaptic transmitter release or a discharge-frequency-calculating interface.
Results
The primary aim of this simulation experiment was to illustrate the advantages offered by this multi-unit model in contrast to the limitations of previous single-unit sphincter models presented in the literature. All simulations were carried out using dedicated software in Borland Delphi Professional version 4 on an Intel Pentium-II computer.
Description of the simulations
The integrated model was used to simulate (i) a normal sphincter and (ii) a partially deficient innervation of the sphincter. The normal situation was simulated with a fully functional neural net providing a uniformly activated sphincter. In the second simulation, one quarter of the sphincter was not activated. Although this activation pattern is unlikely to occur under physiological conditions, it was chosen to test the robustness of the present model. In the experiments, the pressure in the lumen of the sphincter was prescribed to increase sinusoidally from 5 to 60 kPa. All neurons in the network received a current input of 2.5 nA in pulses of 2 ms, with an interspike interval of 55 ms, causing the neurons to fire with a frequency of approximately 20 Hz.
The connection between some neurons and the accompanying muscle elements was blocked to simulate a partially denervated sphincter. Muscle properties did not change between simulations. The results of the simulation are shown in Fig. 3. Fig. 3A–D shows the results of the partial activation experiment at various times (t=0, 350, 700 and 1050 ms), while Fig. 3E shows the state of the normally activated model at 1050 ms.
The network activity is plotted on the top row in Fig. 3. Each rectangle represents a neuron, forming a 3×7 matrix of neurons. Each muscle element is innervated by one neuron. The model sphincter consists of five circular layers. The three outer layers represent muscle tissue, while the inside of the sphincter is formed by two circular layers of passive tissue to simulate the epithelium and the lamina propria (visco-elastic connective tissue). The vertical columns of neurons innervate the circular layers of the muscle, starting from the outside, whereas the horizontal rows innervate the radial segments of the muscle, starting at the 09:00 h position and proceeding in a clockwise direction (see Fig. 2C). This configuration is artificial but can easily be replaced if detailed anatomical data become available. The two right-hand columns of the neural network in Fig. 3A–D show a (prescribed) normalised activity of zero in contrast to the network activity in the normal situation (see Fig. 3E).
The network activity stimulated the sphincter, resulting in an active state of the muscle elements (second row in Fig. 3). A fully activated network resulted in a uniformly activated sphincter (see Fig. 3E), whereas a partially inactivated network resulted in a non-uniformly activated sphincter (see Fig. 3A–D).
The fibre strain (third row in Fig. 3) is low when t=0 ms. The sphincter is pushed outwards by the increased luminal pressure. This causes the strain in the sphincter wall to increase. However, since a portion of the sphincter is activated and tends to contract, only the non-activated portion of the sphincter is stretched. This is clearly shown in Fig. 3D, in which the strain in the activated portion is approximately zero, but the stress is high. Thus, the strain graphs show that the non-activated portion of the sphincter is stretched by the activated contracting portion of the muscle elements and the increasing pressure in the lumen. Counteracting forces act on the sphincter, causing high pressure near the lumen and low pressure on the outside of the muscle. The gradient is proportional to the tensile stress (fifth row in Fig. 3) and is inversely proportional to the radius of curvature of the layers of muscle fibres. The smaller the curvature (at positions nearer the middle point of the sphincter), the higher the pressure gradient that can be generated.
A few remarks can be made on comparing the model states at t=1050 ms in the experiment (Fig. 3D) and the control (Fig. 3E). First, there is no difference in muscle properties; both states reflect the effects of distinct activation patterns on normally functioning sphincter tissue. Second, the strain in the non-activated portion is so high that the non-activated portion bulges out. By changing the uniformity of activation of the muscle fibres, the sphincter loses its normal function.
In conclusion, when the muscle is intact and functioning, the uniformly activated muscle is able to keep the sphincter closed, although the internal luminal pressure is high. Conversely, the partially denervated (non-uniformly activated) muscle is not capable of remaining closed.
Discussion
The integrated neuromuscular model is divided into two sub-models: a neural network and a sphincter muscle.
Neural network
We used a biophysically minimal motor neuron model to simulate neuromuscular sphincter control. In future work, the motor neuron model needs to be extended since Feirabend et al. (1997) have demonstrated the existence of gap junctions in human Onuf’s nucleus, which is the main motor neuron population innervating the urethral and anal sphincters. Gap junctions are of particular importance in forming low-resistance pathways through which currents generated by action potentials in one cell activate adjacent cells. The model offers the possibility of adding gap junctions between compartments of separate neurons. Although the physiological experimental data of Sasaki (1991) were partially implemented, the experimentally indicated anomalous inward rectifying current in external urethral sphincter motor neurons analogous to the hyperpolarization activated IQt current (Halliwell and Adams, 1982) (also known as Ih) (Mayer and Westbrook, 1983) has not been incorporated.
The inferior olive nucleus shows similar characteristics: the cells are coupled by electrotonic gap junctions located on the dendrites (Llinás et al., 1974; Sotelo et al., 1974). Ih was also found in the inferior olive nucleus. Detailed studies on Ih have revealed its critical role in synchronised oscillations (McCormick and Pape 1990; Soltesz et al., 1991; McCormick and Huguenard, 1992; Bal and McCormick, 1997). Consequently, future modelling needs to include these items.
Sensory feedback was implemented to represent a basic reflex system. Nevertheless, during the simulations, the model functions as an open-loop system because the input of the strain receptor activity (=sensory feedback) into the network was blocked to ensure uniform neural activity in the neural network during the simulation experiment.
The simplification of the motor neurons was motivated (i) by the minimal experimental data available about external urethral sphincter motor neurons from Onuf’s nucleus and (ii) by the fact that a more complicated model would make no difference in the present simulation experiment. In this simulation experiment, the motor neurons are responsible only for the uniform and synchronous activation of the sphincter muscle.
Sphincter muscle
The description of a multi-unit muscle as a two-dimensional system of point masses has also been used in a simulation of a non-circular muscle: the tentacle of squid (Van Leeuwen and Kier, 1997). The simulation of the completely activated sphincter model showed a circular muscle. Various references in the literature indeed show circular (ultrasound) images of the morphology of the (anal) sphincter (Law et al., 1990; Tjandra et al., 1992; Gantke et al., 1993; Schäfer et al., 1994). The model presented here is not limited to axi-symmetric problems, although the initial geometry chosen for the modelled structure is. The model is able to simulate a variety of (sphincter) muscle forms, such as a U-shaped geometry. In our model, muscle fibre activation is not mediated by biophysical correlates such as presynaptic transmitter release, postsynaptic ligand-activated conductance, end-plate potential and Ca2+ action potential in T-tubules, but only by the discharge frequency of the motor neuron: activation of a motor neuron results in contraction of the fibres it innervates. Relaxation or lengthening of muscle elements can be achieved passively (by the lumen pressure), for which it is necessary to reduce the motor neuron discharge (Berne and Levy, 1994). The tension generated in a muscle motor unit will, therefore, be a function of the discharge rate of its motor neuron.
The model presented here can handle more advanced realistic sphincter morphologies. The subdivision in separate circumferential muscular elements opens the perspective for realistic representations of, for example, the distribution of activation. Local variations in the sphincter morphology, such as muscle fibre types (Light et al., 1997; Krier et al., 1988; Oelrich, 1983; Schrøder and Reske Nielsen, 1983; Oelrich, 1980; Bowen et al., 1976), muscle and connective tissue fibre arrangements and neural activation can be simulated using the present model. This model cannot deal with muscle fibres traversing the defined layers. Future extensions of the model might support three-dimensional architectures.
The simulation
The experiments described in this paper were simulated further than t=1050 ms, but no significant changes could be detected in the remainder of the simulation. Test runs confirmed that the present model is stable for at least 20 s.
Physical stability
In the simulation of uniform muscle activation, the models evolve to an equilibrium configuration in which all muscle fibres are shortened. The simulation of partial activation is an extreme test of numerical and physical stability. Partial activation leads to a physical state that is initially far from equilibrium. Neighbouring fibres (in a circumferential arc) display completely different behaviour. The activated fibres contract, resulting in stretching of the inactive fibres. The abrupt transition between the activated and the non-activated portion results in a large deviation in axial symmetry and in abrupt jumps in the local strain. The model provides a numerically stable solution of the highly transient behaviour towards the new equilibrium situation (with the applied steady external load). We expect the model to behave equally well in less extreme physical circumstances.
Future verification of the model
At present, detailed experimental data for the strain, stress and pressure distribution in sphincter muscles are rare, presumably because of experimental difficulties. Future research requires advanced experimental techniques to collect the data necessary to verify the results of the model. Promising experimental techniques seem to be high-resolution nuclear magnetic resonance tagging, ultrasound imaging, the application of miniature ultrasound crystals for distance measurements and distributed needle electromyography. The limiting factor of current imaging research is the low resolution; much higher resolution is needed for our purposes. The present model forms the best available prediction of sphincter deformations and may be a useful tool in the choice of the experimental arrangement required to study the spatial behaviour of sphincters. The need for detailed experimental data could direct future experimental research in different areas such as those mentioned above.
However, using advanced models in the absence of experimental data is the second-best option for studying complex non-linear systems with distributed properties. The model presented here points to an interesting phenomenon: partial activation of a sphincter muscle leads to work done on the inactive portion by the activated portion and, consequently, to a less effective sphincter muscle. This simulation illustrates an application and the potency of integrated multi-unit models in contrast to single-unit muscle models.
Discussion of the dynamic stability and the consequences of non-linearity to the control of movement, although beyond the scope of the present paper, is a logical next step. In future research, this type of model supported by detailed experimental data might be helpful in biological and/or clinical applications.
Acknowledgements
The authors are indebted to Professor Dr A. R. Bakker, Professor Dr E. Marani and Professor Dr A. A. B. Lycklama à Nijeholt for useful suggestions.