In vitro virtual reality: an anatomically explicit musculoskeletal simulation powered by in vitro muscle using closed-loop tissue–software interaction

Skeletal muscle is astonishingly versatile. Despite its conserved contractile machinery across taxa (Lindstedt et al., 1998), a single muscle may shift from motor to brake, depending upon the relative timing of force and length patterns (Full et al., 1998; Ahn and Full, 2002; Ahn et al., 2006). Between these mechanical extremes, animals fine-tune muscle force–length dynamics to modulate mechanical power for running (e.g. Daley and Biewener, 2003; Gabaldón et al., 2004), swimming (e.g. Rome et al., 1993; Altringham and Johnston, 1990; Marsh et al., 1992; Richards and Biewener, 2007) or flying (e.g. Askew and Marsh, 2001; Ellerby and Askew, 2007; Morris and Askew, 2010). Our understanding of in vivo muscle versatility stems from in vitro experimentation where isolated tissue is stimulated at various fixed lengths (Ramsey and Street, 1940; Gordon et al., 1966), fixed loads (Hill, 1938) or fixed shortening–lengthening cycles (i.e. work loops; Pringle and Tregear, 1969; Josephson, 1985).

The above approaches are limited because mechanical loading is predetermined. Although changes in muscle stimulation can alter muscle work output, the fixed loading conditions do not permit variation in the acceleration trajectory of the load because length (or force) input is controlled at each time point. Hence, in vitro experiments can ‘replay’ mechanical events to mimic in vivo observations, but cannot explore how hypothetical changes in the neuro-musculoskeletal system might influence movements of the load. For example, rocket frogs achieve extreme jumps via both specialisations in contractile properties and increased tibia length (James and Wilson, 2008). Yet, although controlled-length experiments reveal how muscle properties confer power (Lutz and Rome, 1994), they cannot determine how bone morphology influences performance; changes in strain from concomitant shifts in leg posture, kinematics and muscle loading would be unknown and difficult to predict a priori.

Towards overcoming the predictive limitations of in vitro experimentation, pioneering studies used closed-loop force-feedback electronics to mimic the effects of varying inertial loads in vitro (Lin and Rymer, 1998, 2000) as well as to manipulate neuronal feedback in vivo (Sponberg et al., 2011). More recently, virtual reality environments have been developed into ‘neuromorphic models’ to interact with cadaver tissue (Niu et al., 2017) as well as simplified mechanical models to test effects of inertia, skeletal gearing and elasticity on in vitro muscle (Robertson and Sawicki, 2015; Sawicki et al., 2015). Building on those important ‘virtual loading’ techniques, as well as recent advances in muscle–hardware feedback (Richards, 2011; Richards and Clemente, 2012; Eberhard and Richards, 2018), we introduce in vitro virtual reality (in vitro-VR) as a hybrid of in vitro physiology (e.g. Lutz and Rome, 1994) and musculoskeletal simulation (Delp et al., 2007). In vitro-VR uses closed-loop feedback between simulation software and physiology hardware to cause an in vitro muscle to ‘feel’ immersed in a virtual reality environment. A crucial difference between in vitro-VR and previous methods is the incorporation of multibody dynamics and realistic ground contact forces. Thus, in vitro-VR has the potential to capture anatomical complexity (as does OpenSim) whilst additionally avoiding the pitfalls and oversimplifications of traditional muscle models (see Millard et al., 2019). For demonstration, we tested the effects of muscle stimulation versus altering limb anatomy by modifying the muscle's anatomical origin. We then describe the methodology and limitations as well as present data demonstrating in vitro-VR's utility for recording force–length dynamics in direct response to manipulating the virtual skeleton (e.g. muscle moment arms, limb inertia) towards a better understanding of how evolutionary transformations in limb architecture influence the mechanical demands of muscle.

In vitro-VR is a hybrid of two components, one part software and the other hardware. (A) A 2D musculoskeletal model of a jumping frog is implemented in a software physics engine (MuJoCo; Todorov et al., 2012), which solves the forward dynamics equations to compute joint kinematics in response to joint torque inputs (Fig. 1A). (B) An in vitro muscle–tendon unit (MTU) is mounted to a custom ergometer (Richards and Clemente, 2012) instrumented with a 305C-LR muscle lever system (Aurora Scientific Ltd, Aurora, ON, Canada) which measures force at a controlled MTU length (Fig. 1B). In real time, components A and B form a closed-loop system; force recorded from B is transmitted to A, and converted to a joint torque via a moment arm r (see below), which causes joint acceleration. At each time step, acceleration is integrated to update instantaneous joint angles and MTU length. To close the loop, MTU length is transmitted back to B to update the instantaneous length of the real muscle tissue. Hence, the system uses a single experimental input (electrical stimulation to the muscle) and two primary measurements: muscle force (measured in B) and length change (computed in A).

In vitro-VR involves a chain of intercommunicating digital devices, each with latencies that ultimately limit the update rate of the closed loop. Two main latencies occur: simulation latency, which is the time required to perform forward dynamics calculations, and communication latency, which is the time required to read/write/convert digital information (see below). MuJoCo has been designed to run faster than real time (i.e. simulation latency is less than simulation time step) for robot control, which requires control signals to be computed/optimised predictively (Erez et al., 2013). Owing to the impressive performance of MuJoCo, simulation latency is tiny compared with communication latency, and hence the musculoskeletal simulation remains idle for the majority of the overall loop cycle. As explained below, communication lag is the rate-limiting step; therefore, the timing of digital communication was used to clock and synchronise all processes in the loop.

Our system architecture and software control minimises latencies among the system components, allowing the hardware–software loop to be updated at 1 kHz (see below). For our model, this update rate is sufficiently fast for muscle to smoothly ‘feel’ the resistance of the simulated limb via the amount of length change fed back through the system. For example, if the simulation encountered high resistance (e.g. a limb pulling a massive object), the resulting length change would be nearly zero at each time step. In contrast, for a low resistance (e.g. a limb moving through air), the muscle length change would be large. Importantly, we discuss (below) cases where 1 kHz is not sufficiently fast and suggest how to increase update rate.

In the current implementation, the jumper has a single leg with four planar hinge joints: tarso-metatarsal (TMT), ankle, knee and hip (Fig. 1A). Simplifications were made to allow the jumper to operate with a single extensor muscle (see below). The body segment masses and moments of inertia were determined from micro-computed tomography (µCT) image analysis of an 11.4 g Kassina maculata (walking/jumping frog) specimen as done previously (Porro et al., 2017; Richards et al., 2018). To extend the ankle and TMT joints, the plantaris muscle was modelled as a rigid cable originating at the knee and inserting on the proximal tarsals (Fig. 1A). Importantly, this ‘cable’ is not meant to represent rigid tissue properties, but rather serves as a placeholder which specifies the line of action along which muscle force is applied in simulation. Solid mass-less spheres on the joint centres allow the muscle path to curve smoothly around the joints without penetrating leg segments. At each simulation time step, MuJoCo's kinematics solver computes moment arms at the TMT and ankle. Consequently, the user does not assign moment arms directly; instead, they are computed based on the model's geometry and instantaneous posture. Note, however, that in the current model, the muscle path is constrained to wrap around the surface of a sphere, thus keeping moment arms nearly constant throughout the jumps. Future implementations could explore how time-varying moment arms might influence jumping (Roberts and Marsh, 2003). The simulated plantaris path represents a numerical placeholder for the in vitro muscle tissue. Although the plantaris contains fibres originating proximal to the knee (Dunlap, 1960), for simplicity, we defined the origin point precisely at the knee centre of rotation (i.e. the muscle produces no torque at the knee; we define this as the nominal condition). When the simulation runs, instantaneous in vitro muscle force enters the simulation as a scalar value which is applied along the spatial path of the modelled muscle to produce both ankle and TMT torque.

Although our current model is implemented based on K. maculata, we used muscle isolated from similar-sized Xenopus laevis to develop and demonstrate the in vitro-VR technique. Future use of this technique to evaluate jumping performance in K. maculata will employ muscle isolated from K. maculata specimens. We stress that the current study used a muscle from this non-specialised jumper only to demonstrate the technique. For the present purposes, we wish to illustrate how real muscle (as opposed to a Hill-type muscle model) can be coupled to a simulation to yield interesting biological insights. However, we also note that anuran muscle properties have been shown to vary across species (Astley, 2016); thus, inter-species effects are likely to be important and should be explored in the future.

An adult male X. laevis was weighed then killed using 2 g l⁻¹ solution of tricaine methanesulfonate (MS-222; Sigma-Aldrich, Dorset, UK) followed by removal of the heart, according to a UK Home Office Licence (70/8242). The plantaris longus muscle was removed, mounted to the ergometer and bathed with standard oxygenated saline solution maintained at 22°C and augmented with CaCl₂ (1.3 mmol l⁻¹ Ca²⁺) and glucose (10 mmol l⁻¹). Tissue preparation is detailed in prior work (Richards and Clemente, 2012, 2013). Briefly, using repeated twitch contractions, stimulation voltage was increased until maximum isometric force was achieved. Subsequently, optimal length (24.08 mm) was found by adjusting passive tension until the muscle produced the highest force (approximating the plateau of the force–length relationship; Richards and Clemente, 2013). Prior to experimental trials, the mass of the current frog specimen (30.3 g) was used to scale the mass/inertia parameters of the MuJoCo frog model (11.4 g) using isometric scaling such that the mass of the model matched that of the real frog specimen whilst maintaining proportions. For frog jumps, the muscle was stimulated at maximal voltage at a spike frequency known to produce maximum force for X. laevis (250 Hz; Richards, 2011). Foot contact time for K. maculata ranges between ∼80 and 200 ms (Porro et al., 2017); thus, a stimulation duration of 80 ms was selected based on the lower bound of this observed range. Note, however, that no attempt was made in the current study to maximise muscle force, which would probably occur at longer stimulation durations. The muscle was allowed to rest for at least 5 min between contractions.

In principle, there are many ways to achieve closed-loop hardware–software feedback for in vitro-VR. For the current implementation, two loops were used: a simulation loop (inner loop) running an instance of MuJoCo and a timing loop (outer loop) generating 1 kHz clock signals (Fig. 1C). A schematic diagram showing hardware parts and wiring is shown in Fig. S1. The code for the current implementation can be found from an online repository (https://github.com/frogtronics/InvitroVR).

The simulation loop has the following sequence of events, where MuJoCo performs the following: (1) joint angles are updated (based on integrated accelerations from the previous time step); (2) moment arms and muscle lengths are computed; (3) in vitro force is substituted into Eqns 5 and 6, and hip/knee servo torques are computed; (4) internal forces are calculated; (5) external forces are solved (e.g. ground reaction force); (6) Eqn 1 is integrated to advance the simulation by one timestep (dt=1 ms) to give the current simulation time (t=t+dt); and (7) simulation state (positions, angles) is stored in a buffer for later saving.

When a clock signal is received from the timing loop, the simulation loop returns to step 1 to repeat the computation sequence. In the current application, steps 1–7 are completed rapidly (<<1 ms) before the simulation loop waits idle for a clock signal to return to step 1.

In addition to fast computation, in vitro-VR requires fast communication between analog (muscle ergometer) and digital hardware (PC). To avoid confusion, input/output will be discussed with respect to the musculoskeletal simulation (in the PC), which receives muscle force input and produces a muscle length output. Analog/digital (A/D:D/A) conversion, device communication, synchronisation and timing were achieved with a high performance Arduino (Arduino Due; 84 MHz 32-bit ARM core microcontroller; Arduino, Somerville, MA, USA). A National Instruments NI-6289 data acquisition system (DAQ; National Instruments, Austin, TX, USA) generates the stimulation pulses to the muscle whilst synchronously recording analog data output from the ergometer at the same rate as the loop (1 kHz). The DAQ runs ‘outside’ the control loop; it operates deterministically in parallel with the microcontroller, but is independent from any computation/communication latencies in the loop. Thus, the DAQ could be used as an independent measure of ‘real world’ time. For best performance, MuJoCo should run on a dedicated PC (or processor core) whilst the DAQ and Arduino software run on a separate PC. Additionally, a field-programmable gate array (FPGA; Compact Rio 9074; National Instruments) runs at 10 kHz to perform three tasks. (1) It applies a second-order low-pass Butterworth filter set at 60 Hz and 100 Hz for length and force signals, respectively. (2) The FPGA also performs voltage matching between devices to allow the signals from the Arduino to drive the muscle lever (see Fig. S1). (3) The FPGA checks the validity of the signals and enables/disables the ergometer accordingly (see code repository; https://github.com/frogtronics/InvitroVR).

Most importantly, transmissions among the devices must be synchronised such that the time elapsed in the simulation (t_sim) matches time elapsed in the real world (t_real). Precisely, for every millisecond elapsed in the real world, the simulation must advance exactly one time step, regardless of the computation time of the simulator. Given that loop rate is limited by communication latencies, synchronisation can be achieved by precise timing of data transmission from the microcontroller to the PC. Specifically, the simulation only advances after receiving a valid signal from the microcontroller which is precisely clocked by a built-in interrupt timer. All digital data were packaged as ASCII characters in binary form and tagged with an ‘end-of-transmission’ code, which validates the data then signals the simulation to proceed. Because MuJoCo is not run deterministically in the current implementation, simulation latency can vary depending on processes that are currently running in the operating system. Before conducting experimental trials, loop delay should be measured (Fig. 1D). Most importantly, a lookback test was performed; the simulation transmitted the force signal as output (instead of length) to verify that an analog signal can travel entirely around the loop without being altered by simulation/communication latencies or by conversion errors (Fig. 1E). Additionally, each trial recorded values of t_sim (from MuJoCo) and t_real (ground truth time from the operating system) as separate data streams to verify post hoc that t_sim=t_real. Together, these validations confirm correct timing (i.e. hardware/software running at the same rate), correct synchronisation (i.e. same start time) and fidelity (i.e. no numerical errors in ASCII decoding or transmission beyond precision/resolution limits).

(1) First, the user initialises the system: the DAQ system is enabled, an instance of MuJoCo is started, the microcontroller establishes an RS-232 link to the PC (and vice versa) and the microcontroller initialises a 1 kHz clock. Once the system is enabled, the microcontroller begins listening for analog force data from the ergometer as well as for valid digital length data from the PC. The simulation is paused at t_sim=t_real=0. (2) The user presses a hardware trigger to initiate the timing loop. (3) The microcontroller converts an analog force signal (from the ergometer force output) to digital then writes the data to the serial RS-232 port of the PC. (4) The PC continuously reads data from the incoming RS-232 port until valid force data are received from the microcontroller then decoded from binary ASCII characters to a floating point force value. This event triggers a single iteration of the simulation loop. (5) The simulation loop advances (see sequence above). (6) The PC converts the instantaneous MTU length (computed from the simulation) to digital and writes the data to the outgoing RS-232 port. (7) The microcontroller reads the MTU length data and decodes it to a floating point value then scales the value using a calibration factor and offset appropriate for the ergometer hardware. (8) The microcontroller writes the MTU length value to the analog output, which is transmitted to the length input of the ergometer hardware. (9) If t_sim is less than the specified end of simulation (e.g. 0.5 s), the loop returns to step 3. (10) After completion of the simulation, a reset sequence is initiated which automatically returns muscle length to its initial value and resets all software values to their nominal values. (11) Buffered simulation data are saved to file. The DAQ saves all data recorded from analog channels.

To demonstrate the method, we performed five arbitrary experimental conditions: (1) a nominal condition where the muscle was stimulated 20 ms following the onset of hip and knee extension, (2) a control condition (control 1) where the plantaris muscle was absent from the simulation, (3) a passive condition (control 2) where the muscle was included in the simulation, but not stimulated (and thus could produce passive force), (4) a nominal condition but with stimulation onset beginning 60 ms following proximal joint motion and (5) a modified anatomical condition where the virtual origin of the plantaris (in simulation) was switched to a proximal location on the femur to cause knee flexor action from plantaris contraction. The nominal condition is based on a well-established observation that frogs extend their proximal joints earlier in jumps (e.g. Calow and Alexander, 1973). Control 1 and 2 are theoretical conditions to verify that the simulation is dependent upon feedback from the real muscle. The conditions for delayed stimulation and modified anatomy (conditions 4 and 5, above) are arbitrary hypothetical conditions to demonstrate how the technique is sensitive to changes in neural stimulation as well as changes in anatomy. Across all conditions, the stimulation burst and the hip/knee kinematics were held constant.

A fundamental limitation of in vitro muscle measurements is that traditional stimulation of the whole muscle or whole nerve cannot replicate in vivo activation where motor units are independently activated according to the mechanical demands of a task (Hodson-Tole and Wakeling, 2009). In the future, one could incorporate sophisticated stimulation electrode designs (e.g. Branner et al., 2001) towards a more physiological stimulation pattern. All other limitations are with regard to anatomical complexity. Our system design only allows a single muscle to be tested in a given closed-loop system. However, additional loops could be easily linked in parallel (each with its own dedicated processor core). This would allow multiple in vitro muscles (each mounted to its own ergometer) to dynamically interact as they each attach (virtually) to a simulated limb. An additional limitation is the lack of tendon dynamics. In future implementations, the user can attach in-series tendons to in vitro muscles by adding spring elements emulated in software. Finally, the present implementation assumes that muscle fascicle strain is identical to whole-muscle strain. In the future, one could use sonomicrometry to measure and control fascicle length using an inner control loop (Robertson et al., 2017).

This section discusses the principal system limitations of the current implementation and suggests modifications for improved performance for future implementations of in vitro-VR. The four principal limitations to in vitro-VR are: (1) maximum speed of the analog ergometer hardware, (2) minimum/maximum numerical range and resolution (i.e. digits of precision for floating point force/length values), (3) loop rate (1/maximum latency) and (4) simulation rate. Firstly, the ergometer motor (Aurora 305C-LR), depending on damping, has a maximum speed of ∼400 mm s⁻¹, which far exceeds the maximum speed of muscle contraction in most applications, and thus is a negligible limitation. For example, a 2 cm muscle will reach 200 mm s⁻¹ at a contractile speed of 10 lengths s⁻¹. In rare cases where faster lever movements are required, the system can be tuned with higher gain settings (at the expense of damping). Secondly, both numerical precision and minimum/maximum range are limited by the resolution of A/D:D/A converters in the loop. Currently, we use the built-in converters on the microcontroller, which operate between 0 and ∼2.2 V (after accounting for a DC offset). The converters are 12-bit allowing the range of 2.2 V to be divided into 4096 levels, giving a numerical precision of 2.2/4096 V, which is sufficiently fine for the current application. Given that our ergometer was calibrated at 1 V=1 N, our muscle force limit was 2.2 N, which was sufficient for the size of the muscle used currently. However, for applications with larger muscles producing force in excess of 2.2 N, either the force would require attenuation or an alternative microcontroller would be used. Thirdly, as stated above, communication latency dominates our implementation. At each time step, the microcontroller/PC must use serial communication to transmit/receive 12-bit values. The speed bottleneck of our current implementation is the PC operating system (64-bit Windows 7) which permits a maximum RS-232 baud rate of 115,200 bits s⁻¹. Combining all delays (four read–write operations, MuJoCo computations and other operations), the maximum overall loop rate of the current system is between 1 and 2 kHz; thus, we used 1 kHz for safe operation. If the system were too slow and/or of low resolution, the muscle would ‘feel’ coarse and grainy motion versus smooth. Additionally, force would oscillate, which can be seen in slower feedback systems (see fig. 2A from Richards, 2011, which suffered from a 2 ms delay). Neither ‘graininess’ nor oscillations were observed in the current implementation, indicating that 1 kHz is sufficiently fast for loads produced at the time scales of a frog jump. Finally, simulation rate must be considered. Although the current implementation is relatively simple (four joints and one contact), more complex simulations could increase simulation latency which could invalidate the dynamics because muscle force would not correspond to simulation acceleration. In principle, additional degrees of freedom increase the complexity of the multibody dynamics. In practice, MuJoCo can handle extreme anatomical complexity (at least 1000-fold greater complexity than the current model) whilst maintaining a simulation latency less than 1 ms (see http://www.mujoco.org/performance.html). We note, however, that contacts and collisions occurring at rapid time scales may significantly slow computation time. For example, we speculate that a large quadruped in a gallop could be extremely challenging to simulate with sufficiently low latency. Future modelling studies would be required to evaluate how model and contact complexity influences simulation latency.

In future implementations, the system performance can be improved. Although 1 kHz allows smooth operation for the current model, it would probably be too slow to model tendon recoil dynamics or for faster dynamics of extremely small frogs. The rate could be increased easily with faster serial communication afforded by increasing the baud rate or by using an alternative communication protocol (e.g. ethernet; Eberhard and Richards, 2018). A faster processor and use of multiple processor cores would increase the speed by orders of magnitude (see http://www.mujoco.org/performance.html). If resolution were limiting (not currently the case), higher quality converters (e.g. 16-bit or higher) would massively increase the numerical precision along with better electronics to improve signal to noise ratio. For future implementations, a deterministic real-time operating system (RTOS) would provide a crucial improvement of performance. Otherwise, the user cannot precisely control the priority of computational processes; an unrelated process (e.g. a software update or virus scan, etc.) may slow computations or communication. Because a true RTOS was not used for the present study, we were required to verify that the control loop was timed correctly to maintain a constant loop rate (Fig. 1D,E). If other processes were running, the loop rate would potentially be time varying rather than constant, rendering the simulation results difficult to interpret.

We measured force and length changes of a muscle responding to simulated inertial and ground reaction forces of a jumping frog within a real-time simulator. Data for all trials are given in Table S1.

Simulated jump trajectory was highly dependent on muscle stimulation conditions as well as muscle anatomy (Fig. 2). To our knowledge, this is the first study to connect an in vitro muscle to an anatomically explicit musculoskeletal simulation using closed-loop feedback. For simplicity and to introduce the technique, we used a single muscle with a 2D simulation. No attempts were made to adjust experimental conditions to produce maximal jumping performance. Future implementations can easily incorporate 3D dynamics and additional muscles ‘in the loop’ (see Materials and Methods) as well as controllers to optimise joint kinematics.

At the onset of all jumps, the simulation initially sank towards the ground as a result of gravity causing a brief period of joint flexion and muscle lengthening (Fig. 3). Soon thereafter, hip, knee, ankle and TMT joint accelerations acted to extend the leg to propel the body upwards and forwards as the plantaris shortened. Two controls were used to verify that real-time muscle simulation feedback influenced jump performance. In control 1 (muscle absent) the leg extended, but jumped only weakly (Movie 1; Fig. 2C). In this case, intersegmental forces among the body segments caused the ankle to extend in the absence of muscle force. This transfer of motion among segments is a crucial, but often under-appreciated phenomenon of multi-body dynamic systems (Kuo, 2001) which would be missed if the mass/inertia of the leg were neglected using simpler modelling and in vitro approaches. Motion transfer occurs as a result of intersegmental forces at the joints; these forces implicitly act as constraint forces attaching each segment to the next (i.e. allowing rotation, but preventing translation). In control 2 (muscle present, but unstimulated), initial joint flexion caused muscle stretching, which produced passive force and ankle extensor torque. This passive force acted like an elastic cable to transmit intersegmental forces from proximal joint extension to the ground. Consequently, the foot produced sufficient ground reaction force to lift the frog farther from the simulated floor compared with control 1 (Movie 2). For subsequent trials with electrically stimulated muscle, the simulation jumped at least twice the height observed for control jumps (Figs 2D and 3; Movie 3). These trials confirm that the closed-loop real-time feedback allows the muscle to influence the simulation and vice versa.

To further demonstrate the utility of in vitro-VR, we altered musculoskeletal anatomy by altering the origin of the plantaris muscle (Fig. 2B). In the nominal case, the plantaris origin was positioned at the knee to eliminate any action about this joint (see Materials and Methods). We explored the effects of editing musculoskeletal architecture by shifting the origin proximally to a point above the knee. This manipulation allowed the plantaris to act at the knee as a flexor, causing dramatic alterations to plantaris force–length dynamics. Rather than a brief period of lengthening followed by a large degree of shortening, the muscle remained nearly isometric then shortened shallowly (Fig. 3B). As the knee extended it acted as an antagonist to the plantaris such that the shortening that would have occurred was instead cancelled by the upward rotation of the femur to ‘pull’ on the muscle origin. Consequently, the plantaris briefly acted as a strut to transmit force to the ankle. In this special case, knee extension is ‘transferred’ to the ankle via a bi-articular link, allowing the ankle to extend with minimal plantaris shortening. Exploring the full dynamic consequences of this shift in force–length dynamics is beyond the scope of the current work as it requires further investigation using a range of incremental anatomical changes as well as a sweep of stimulation conditions. Regardless, for the purposes of demonstrating the technique, in vitro-VR offers intriguing direct evidence that changes in musculoskeletal architecture can greatly influence muscle function, independently of neural stimulation.

In vitro experimentation, particularly the work-loop method (Pringle and Tregear, 1969; Josephson, 1985), enables exploration of parameters to map muscle ‘performance space’ relating force, work and power to strain and stimulation patterns (Josephson, 1985; Altringham and Johnston, 1990; Luiker and Stevens, 1993). One can then ask further questions: why do muscles operate within a particular region in their performance space?; why do muscles perform maximally under some conditions (e.g. Rome et al., 1993), but sub-maximally under others (Tu and Daniel, 2004)? These questions require modelling the physical forces resisting muscle contraction (Marsh, 1999). In some cases, simple modelling can demonstrate how ‘environmental feedback’ (Aerts and Nauwelaerts, 2009) can reduce the ‘volume’ of muscle performance space (Richards, 2011; Richards and Clemente, 2013; Clemente and Richards, 2013). In other cases, for example in swimming fish, sophisticated modelling is required to reveal how optimal muscle function is dependent upon the mechanical properties of the body, which governs fluid–structure interactions (Tytell et al., 2010). Unfortunately, however, modelling struggles to capture the full richness of muscle function observed during active lengthening (as reviewed by Herzog, 2014; Nishikawa, 2016) or extremely rapid events such as ‘short-range’ length perturbations (Millard et al., 2019). Consequently, traditional Hill-type models will fail to explain rapid force–velocity transients occurring on millisecond time scales crucial in behaviours such as ballistic tongue projection in frogs (Lappin et al., 2006) or the response to unexpected perturbations in running birds (Daley et al., 2006, 2009).

In light of the above shortcomings of muscle modelling approaches, there are large gaps between our understanding of in vitro contractile dynamics and in vivo function at the organismal level. We propose in vitro-VR can help close this gap by measuring muscle force–length dynamics in response to an arbitrarily realistic model of limb anatomy, inertia and substrate contact. In essence, in vitro-VR directly measures aspects that remain poorly characterised and are challenging to model (i.e. tissue dynamics) whilst simulating the remaining parameters that can be modelled and validated more easily (multibody and contact dynamics). We make a crucial point that in vitro-VR does not replace traditional in vitro experimentation; the methods are complementary. In vitro experiments such as work loops reveal the full performance space of muscle whilst in vitro-VR explains how anatomy and environment limit in vivo operation within this space.

In vitro-VR can provide unique insights into several current problems. Firstly, several researchers have refined the Hill-type muscle model (Gordon et al., 1966; Hill, 1938); however, muscle models are difficult to validate outside of standard protocols where parameters are held fixed. In vitro-VR can produce contractile dynamics closer to in vivo conditions because the dynamic loads relate directly to limb function observed in natural behaviour (i.e. activation, velocity and force each vary simultaneously). Thus, in vitro-VR could help validate recent muscle models accounting for traditionally neglected features such as the role of titin to explain history-dependent effects (Nishikawa et al., 2011) as well as ‘short-range’ impedance effects (Millard et al., 2019), tissue mass (Günther et al., 2012; Ross and Wakeling, 2016) and 3D structure (Blemker and Delp, 2005; Ross et al., 2018). Secondly, in vitro techniques have established how force is enhanced through stimulation timing (e.g. Burke et al., 1970; Sandercock and Heckman, 1997), active stretching (Edman et al., 1978) and starting length (Azizi and Roberts, 2010). In vitro-VR can be used to determine whether these well-established effects are likely to occur in vivo and whether they enhance performance during a complex behaviour. Thirdly, in vitro-VR can be an important tool to understand how muscle–tendon systems ‘amplify’ muscle power for propulsion (Peplowski and Marsh, 1997) or ‘attenuate’ power for braking (Ahn et al., 2006; Richards and Sawicki, 2012; Konow and Roberts, 2015), depending on loading and posture (Galantis and Woledge, 2003; Roberts and Marsh, 2003) as well as substrate properties (Reynaga et al., 2019). Our method could further elucidate how dynamic inertial effects of multi-body systems (Richards, 2019) influence these loading and substrate effects. Finally, in vitro-VR could be used to directly measure the functional consequences of anatomical changes occurring over evolutionary time (e.g. Emerson, 1979; Reilly and Jorgensen, 2011) or due to ecological factors (e.g. Louppe et al., 2017).

We present a technique, in vitro-VR, which aims to advance in vitro muscle physiology and musculoskeletal simulation (e.g. Delp et al., 2007) by creating a novel link between the two approaches. In vitro experimentation enables physiologists to control neuromuscular stimulation and strain patterns to help determine how the nervous system modulates mechanical work and power. In addition to stimulation, in vitro-VR allows manipulation of limb anatomy and locomotor substrate. When traditional methods and in vitro-VR experimentation are applied in tandem, one can probe how anatomy and locomotor environment influence the mechanical demands on a muscle in vivo. One can address new questions, such as: how does anatomy or substrate influence the mapping between motor stimulation and mechanical output? Additionally, the present method has wider applications; for example, to test tissue samples affected by disease or injury to better isolate and understand their effects on locomotion. Hence, we propose in vitro-VR is a broadly applicable method within and beyond the comparative/evolutionary context.

We greatly thank Tim West for advice on muscle preparation, and assistance with laboratory setup and buffer solutions. We also thank Rue Jones for donating to us a small number of Xenopus frogs for this study. Cameron Hill assisted with equipment testing and Chris Basu and two anonymous reviewers provided helpful comments on the manuscript.