{'fraction_non_alphanumeric': 0.031342975815241404, 'fraction_numerical': 0.010899343093645977, 'mean_word_length': 4.719966301600674, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neuromuscular actuation of biohybrid motile bots

Significance Biohybrid machines have been developed using muscles to actuate soft robotic structures. We envision that the integration of neurons into the embodiment of such systems can transform them into intelligent machines which, for instance, could use sensory neurons to detect environmental cues then adaptively respond and orchestrate various motor patterns through their neural circuitry. However, achieving such sensory-motor modalities relies on the ability of neural units to command muscle activity, making actuation through motor neurons the first milestone. Here, we achieve this milestone and demonstrate neuromuscular actuation of a biohybrid swimmer. This paves the way for the development of biohybrid embodied platforms as models to gain deeper understanding of motor control, with potentially broad impact in robotics, bioengineering, and health.


Fabrication and assembly of PDMS scaffolds
PDMS structures were cast from microfabricated silicon molds. Silicon wafers were patterned by photolithography, etched using the Bosch process, and coated with polytetrafluoroethylene to facilitate removal of PDMS from the mold. Wafers were etched to a nominal depth of 100 µm for the co-culture test platforms and swimmer scaffolds, and to 300 µm for the spacers and muscle seeding molds. PDMS (Sylgard 184) base and cross-linker were mixed at 10:1 ratio by weight, pipetted into the molds, and allowed to spread. Trenches for slender features such as the legs and tails were filled by capillary mircomolding [3]. Samples were cured at 60 • C for 12hr and peeled off the silicon mold. The parts were then assembled on top of a glass coverslip with the co-culture scaffold placed on the spacers and sandwiched between the two muscle seeding molds. Scaffold assemblies were sterilized by autoclaving.

Tissue seeding on test platform
For all tissue seeding processes, ECM solution was prepared on ice by first neutralizing type I collagen (Corning) with sodium hydroxide, and then mixing neutralized collagen thoroughly with Matrigel (Corning). Collagen and Matrigel were used at final concentrations of 2 mg/ml each. Fig. S1 shows a step-by-step illustration of the tissue seeding procedure. To form muscle strips, C2C12 myoblasts were suspended in ECM solution at a density of 2.5 × 10 6 cells/ml. Cell-ECM mixture was pipetted into the cavity formed by the muscle seeding molds and polymerized at room temperature for 30 min. Samples were then inundated in muscle growth medium and incubated for 2 days. At day 2, the muscle seeding molds were removed, parts of the scaffold were cut out to allow free deflection of the legs, and culture medium was switched to muscle differentiation medium. Samples were kept in muscle differentiation medium for 6 days with daily medium replacements. To seed the neurosphere, culture medium was aspirated, fresh ECM solution was pipetted onto the scaffold, then a neurosphere with approximately 400 µm diameter was selected and seeded into the hollow cavity on the scaffold. ECM solution was then allowed to polymerize at room temperature for 30 min after which the samples were inundated in co-culture medium. Co-culture medium was replaced daily until experiments were terminated.

Swimmer fabrication, tissue seeding, and release
The swimmer scaffold design included several modifications of the test platform, and the tissue seeding procedure was adjusted accordingly (Fig. S4). The head and legs of the swimmer have the same geometry as the test scaffold except for the swimmer head being rounded to reduce drag. Slender tails emanate from each leg and are attached at the ends to the scaffold via compliant springs. Neuron-muscle co-culture was developed on the swimmer scaffold (Fig. S4) following a similar procedure as for the test platform. To achieve untethered locomotion, the swimmer was released by cutting around the head and at the ends of the tails. The swimmer was immersed in 10 ml of co-culture medium in a 100 mm tissue culture dish. Then, 10 ml of a 1:1 mixture of Percoll (Sigma Aldrich) and co-culture medium was carefully injected at the bottom of the dish to suspend the swimmer in fluid. To confirm that the swimmer is suspended to a sufficient height to prevent friction, we took z-stack images and measured the distance between the swimmer and the bottom of the dish. The bottom of the dish was scored manually using a glass pipette tip to identify when the bottom surface comes into focus. Z-stacks were taken at 10 µm intervals starting from below the dish up to a plane above the swimmer (Movie S2). Before adding Percoll-medium mixture, the swimmer was at a height of approximately 250 µm, comparable to the thickness of the swimmer scaffold (100 µm). After adding 10 ml of Percoll-medium mixture, we allowed 30 min for the flow perturbations to diminish and then took z-stacks every 15 min for 6 hr. During this time, the swimmer remained suspended at a height of approximately 2 mm which is 20 times the thickness of the swimmer (Fig. S7).
Image acquisition, processing, and optical stimulation All imaging was performed on an Olympus IX81 inverted microscope (Olympus America) with a digital CMOS camera (Hamamatsu), mounted on a vibration isolation table. The microscope was equipped with an environmental chamber to maintain samples at 37 • C and 5% CO 2 during imaging. To visualize motor neurons, fluorescent images of Hb9-GFP+ MNs were taken using a GFP filter coupled to an X-Cite 120PC Q widefield fluorescent light source (Excelitas Technologies). To record muscle contractions and swimming, brightfield images were taken at 100 fps using a 4X air objective. Leg deflections caused by muscle contractions and the swimmer's position during swimming were measured from video recordings using the image processing software Tracker (http://physlets.org/tracker). Optical stimulation was performed using the fluorescent light source and GFP filter, delivering blue light at 470 nm wavelength with 3.9 mW/mm 2 intensity as measured by a power meter at the sample plane. Samples were stimulated with 200 ms or 1000 ms exposure by controlling the shutter of the fluorescent light source.

Correction of swimmer position for drift
While the swimmer is kept free-floating in a dish, it drifts relatively slowly due to fluid flow caused by manual handling of the petri dish. We account for this drift by correcting the position data that is obtained from videos. We recorded videos of the swimmer continuously during the rest phase and swimming phase. Given that the diameter of the dish (100 mm) is much larger than the body size of the swimmer (∼ 1 mm), we assume that the swimmer experiences a uniform and constant drift flow. Thus, the scalar components of drift velocity u x,drif t and u y,drif t were computed from linear fits to the corresponding position component-time data during the 10 s drift period before initiation of active swimming. The components of position data during swimming were then corrected as where x(t) and y(t) are the original data, x xor (t) and y xor (t) are corrected positions, and t start is the time of initiation of active swimming. The net swimmer position is computed as r cor (t) = x 2 cor (t) + y 2 cor (t) and r(t) = x 2 (t) + y 2 (t) for the corrected and original data, respectively. Swimming velocity is computed from the linear fit to the position-time data during the steady-state swimming phase. Fig. S8 shows a comparison of the drift-corrected and uncorrected data, as well as simulation prediction, for the experiment shown in Fig. 4D in main text.

Modeling and simulation of Cosserat rods
We mathematically describe a slender rod (Fig. S2b) by a centerlinex(s, t) ∈ R 3 and a rotation matrix Q(s, t) = {d 1 ,d 2 ,d 3 } −1 , which leads to a general relation between frames for any vector v: v = Qv,v = Q T v, wherev denotes a vector in the lab frame and v is a vector in the local frame. Here s ∈ [0, L 0 ] is the material coordinate of a rod of rest-length L 0 , L denotes the deformed filament length and t is time. If the rod is unsheared,d 3 points along the centerline tangent ∂ sx =x s whiled 1 andd 2 span the normal-binormal plane. Shearing and extension shift d 3 away fromx s , which can be quantified with the shear vector σ = Q(x s −d 3 ) = Qx s − d 3 in the local frame. The curvature vector κ encodes Q's rotation rate along the material coordinate ∂ s d j = κ × d j , while the angular velocity ω is defined by ∂ t d j = ω × d j . We also define the velocity of the centerlinev = ∂ tx and, in the rest configuration, the bending stiffness matrix B, shearing stiffness matrix S, second area moment of inertia I, cross-sectional area A, and mass per unit length ρ. Then, the dynamics [4] of a soft slender body is described by: where Eqs. S1 and S2 represent linear and angular momentum balance at every cross section, e = |x s | is the local stretching factor, andf and c are the external force and couple line densities, respectively (shown in Fig.S2c). The external force and couple line densities can then be conveniently employed to model a number of biological and physical effects, such as muscle activity, connectivity among multiple rods, contact forces, and hydrodynamic loads, as further described in the following sections.
For numerical implementation, the continuous representation given above is discretized into n + 1 nodes of position x i and n connecting cylindrical segments (Fig.S2d), so that midline linear displacements are determined by the internal and external forces acting at the nodes (of mass m i ), while rotations are accounted for via couples applied to the cylindrical elements. The dynamic behavior of a rod is then computed by integrating the discretized set of equations in time via a second order position Verlet scheme. The details of our numerical implementation can be found in Ref. [4], together with a rigorous validation against a number of benchmark problems with known analytic solutions as well as experimental investigations involving contact, anisotropic surface friction, and highly viscous fluids.

Computational modeling of swimmer scaffold
Microfabrication of the swimmers yields structures with rectangular cross-sections. For numerical compatibility with our simulation approach, these structures need to be converted into rods (Fig. S2a) while conserving their mass and moment of inertia. Thus, the following equations are applied to convert rectangular beams into rods: where b, h and E beam are the width, thickness, and Young's modulus of the rectangular beams, and r and E rod are the equivalent radius and Young's modulus of the rod, respectively. Eq. S3 ensures that cross-sectional area and therefore mass per unit length is retained (material density kept constant). Eq. S4 is then used to compute the equivalent Young's modulus for the rod such that the moment of inertia is conserved. We assign the value of the Young's modulus of PDMS to E beam , and use the density of PDMS (Table S2). To account for the different thicknesses of different parts of the scaffold (i.e. tails, legs, and head), the scaffold is modeled as a rod with varying radii, while retaining the mass and moment of inertia of the original structures.

Computational modeling of muscle
We simulate the muscle tissue as a bundle of Cosserat rods which allows us to use the same theory to model the entire biohybrid system. To incorporate muscle activity into our model, the muscle-generated force is computed and wrapped into the external force vectorf in Eq. S1. We write the total force generated by muscle contractility (F m ) as: where F out is the force that muscle tissue exerts on the scaffold ("muscle output"), characterized by the resulting leg deflections, E m is the Young's modulus of muscle, A 0 is the muscle crosssectional area at rest, and is the axial strain. The second term on the right-hand side of Eq. S5 takes into account the force that is exerted to deform the tissue itself. At any time, this force can be written as a product of the internal elastic stress E m (Hooke's law) and the instantaneous cross-sectional area A. However, because the cross-sectional area is measured at rest, we derive a relation between instantaneous area and area at rest as follows: Since muscle tissue is nearly incompressible [5], one can write a volume conservation relation AL = A 0 L 0 . Using L/L 0 = 1 − (strain is taken positive for contraction), we write A = A 0 /(1 − ), hence the form of the second term on the right-hand side of Eq. S5.
As described in the main text, the force muscle exerts on the scaffold (F out ) consists of quasistatic rest tension and time-dependent active contraction force. We therefore write F out as a sum of these two components: where F rest and F peak are the magnitudes of rest tension and peak active contraction force, f is the muscle contraction frequency, and t is time. We approximate the shape of a single contraction as a Gaussian function e −Ct 2 where C is a constant that adjusts the width (duration) of contraction. To simulate periodic contractions with frequency f , we modify this term to N n=1 e −C(t−(n−0.5)/f ) 2 , which gives a series of N contractions where the time between two subsequent contractions is 1/f . The values of F rest , F peak , and f are experimentally determined to be 17.5 µN, 7.64 µN, and 4.2 Hz, respectively, from measured displacement dynamics of the legs, given the known E and geometry of the test platform. The value of C is determined to be 556 by fitting to experimentally measured dynamics of leg deflections. We have demonstrated that when this muscle model is incorporated into our coupled simulations, the resulting leg deflections match the experimental data very closely (Fig. 3B).

Interactions between rods
Full simulation of the swimmer involves assembling multiple rods. To avoid interpenetration between the rods representing the muscle bundle and the scaffold, we introduce repulsive forces F c acting between rods that are in contact. Considering two rods (rod i, j, with radii r i and r j ) that are in contact, the corresponding F c is stated as: where d ij min is the minimum distance vector between the two rods, so that the scalar overlap displacement is written as ij = (r i +r j −d ij min ). The first term within the square brackets expresses the linear response to the interpenetration distance as modulated by the stiffness k, while the second damping term models contact dissipation and is proportional to the coefficient γ and the interpene- is the Heaviside function that ensures a repulsion force is produced only in case of contact ( ij ≥ 0).

Hydrodynamics
As discussed in the main text, due to the targeted low Re operating regime and slender geometry of the swimmer tails, we apply slender body theory to model the hydrodynamics of the tails. The hydrodynamic force per unit length (f H ) acting on the tail is captured by: where µ is the dynamic viscosity, L and r are the length and radius of the rod representing the tail, I is the identity tensor, and v and t denote the local velocity and the tangent vector along the body, respectively. We confirmed the validity of this approach for modeling the elastohydrodynamics of our tails by comparison to benchmark finite element analysis (discussed in the following section).
Motion of the swimmer is resisted by a viscous drag force, primarily acting on the relatively large assembly of head, legs, and neuromuscular tissue that engulfs the legs and part of the head. We approximate this irregular-shaped body with a solid sphere and compute the drag force (f D ) using Stokes' law where R is the radius of the sphere, µ is the fluid viscosity, and v c is the swimming velocity. Although Eq. S9 is derived for a sphere in low Re flow, drag forces on non-spherical objects can also be accurately predicted by the same equation multiplied by a scalar (often empirically determined) correction factor [6]. This is mathematically equivalent to tuning the value of R. For our swimmer, the in-plane dimensions of the head-legs assembly was 1 × 1.5 mm and the thickness of the engulfing neuromuscular tissue was roughly 0.5 mm. We therefore used a sphere with a diameter of 1 mm as an approximation.

Validation against finite element analysis
To demonstrate the capability of our model to predict the elastohydrodynamic behavior of the swimmer tails, we present here three benchmark studies where we validate our model predictions against finite element analysis (FEA) simulations carried out using ANSYS. FEA simulation settings are provided in Table S2. In all FEA simulations, the tail was modeled as a rectangular beam with 10 × 100 µm cross-section (similar to the fabricated tails). In the corresponding Cosserat rod models, as before, equivalent rod radius and Young's modulus were computed using Eqs. S3 and S4. The FEA beams were assigned the Young's modulus of PDMS (Table S2).
As a first step, we performed an analysis of the frequency response in vacuum for tail lengths of 1 mm and 2 mm. The base of the tail was constrained in a plane perpendicular to the long axis of the tail (y-axis) and allowed to move only in the z-direction (see Fig. S3a). The opposite end of the tail was free. Prescribed sinusoidal displacement was applied at the base along the z-direction, with amplitude 10 µm and frequency varying from 1 to 160 Hz, to cover a range that includes the theoretically predicted natural frequency. The resulting tail vibrations were characterized in terms of the amplitude of oscillation of the free end, and plotted against excitation frequency for the two different tail lengths. The results show very close match between the predictions of our Cosserat rod-based models, FEA results, and theoretical solutions (Fig. S3b). This benchmark study confirms that our modeling approach accurately captures the elastic and resonant response of the tail structure.
Next, we performed two case studies to validate our model predictions in terms of the tail elastohydrodynamics. In both cases, slender body theory was used to compute hydrodynamic forces on the Cosserat rod models (as described above), and FEA was carried out using high fidelity, 2-way coupled fluid-structure interaction (FSI) simulation with the beams immersed in a fluid domain. In the first study, boundary conditions were set to mimic the experimental operation of the tails: The base of the tail was constrained in-plane and excited by applying prescribed displacement in the z-direction (triangle waveform, 20 µm amplitude, 1.25 Hz frequency), while the opposite end of the tail was free (Fig. S3c). As a further validation, we performed a second case study with different boundary conditions: Here, the base of the tail was fixed (zero displacement in all directions) and the opposite end was excited by applying prescribed force (triangle waveform, 1.5 nN amplitude, 1.25 Hz frequency) in the z-direction (Fig. S3f). FEA results for beam deflections and flow velocities around the beam are illustrated in Fig. S3d and Fig. S3g. We evaluated tip displacements of the tails and compared the predictions of FEA simulations and Cosserat rod-based models. For both case studies, results showed very close match, confirming the validity of our modeling approach.
To ensure that the hydrodynamic operating regime in the simulated cases were similar to that of experiment, we computed the "frequency Reynolds number" (Re f ) for each case: For a tail oscillating in a fluid (heaving motion) with oscillation amplitude (displacement) A and frequency f , the characteristic velocity is given by Af , and the hydrodynamic operating regime can be described by a frequency Reynolds number [7,8] defined as Re f = ρAf L/µ, where ρ and µ are the density and dynamic viscosity of the fluid, and L is the tail length. Results (Table S3) show that Re f for all simulation cases were on the same order of magnitude as that of experiment.

Optimizing swimming performance
At low Reynolds number (Re 1), inertial forces are negligible and viscous forces dominate. Thus, the locomotion of a low Re swimmer, driven by the thrust generated by the propulsion mechanism, is resisted predominantly by the viscous drag force. Optimizing swimming performance therefore involves two sets of strategies, which aim at maximizing thrust and minimizing drag, respectively. To explore possible strategies for optimizing thrust, we begin with a physical description of propulsion at low Re: Due to the dominance of viscous forces over inertial forces, mechanisms of propulsion that rely on inertia do not work, and the swimmer must continually generate body deformation to maintain thrust. Furthermore, a swimmer operating at Re 1 experiences a flow regime (Stokes flow) for which the governing equations of motion are timeindependent. An immediate consequence of this, known as Purcell's "scallop theorem" [9], is that a swimmer that generates reciprocal (time-reversible) movements cannot achieve net locomotion at low Re since the thrust forces produced as the swimmer's body goes from rest to the deformed configuration are canceled out by forces generated while the swimmer's body goes back to the rest configuration -regardless of the relative speeds with which the swimmer goes through these cycles. As Purcell pointed out [9], one way to achieve time-irreversible deformation is to actuate a flexible filament (by rotating or translating at one end) in a viscous fluid. Indeed, many microorganisms that swim at low Re employ flexible slender appendages (flagella and cilia) to propel themselves [10,11,12].
In our demonstration, we utilized flexible PDMS filaments attached at their base to the legs of the swimmer scaffold. Cyclic muscle contractions evoked by neuronal stimulation induce leg deflections and thereby actuate the filaments by applying displacement at the base (Movie S1). In general, when a flexible and slender filament in viscous medium is actuated at its base, this leads to the propagation of bending waves along the length of the filament [13,14,15,16] which leads to time-irreversible motion. Together with the drag anisotropy (slender body theory gives: ζ ⊥ ≈ 2ζ where ζ ⊥ and ζ are the viscous drag coefficients perpendicular and parallel to the filament axis, respectively), this time irreversible motion can generate net thrust [17]. Under the assumptions of low Re, slender body, planar motions, and small displacements (displacements perpendicular to the filament filament length), and for a given actuator design, the elastohydrodynamics of a flexible slender filament is characterized by the nondimensional Sperm number, where L, E, and I are the length, Young's modulus, and moment of inertia of the filament, respectively, ω is the actuation frequency in rad/s (ω = 2πf ), and ζ ⊥ is the viscous drag coefficient perpendicular to the filament axis. ζ ⊥ is given in slender body theory [18] by ζ ⊥ = 4πµ/[ln(L/r)] where µ is the dynamic viscosity and r is the filament radius (or equivalent radius for a non-circular cross-section).
A physical interpretation of the Sperm number is as a ratio of the drag force acting perpendicular to the filament as it bends (∼ ζ ⊥ ωL per unit length), to the elastic restoring force ("bending force") which acts to straighten the filament (∼ EI/L 3 per unit length). In the limit Sp 1, bending forces dominate, the filament moves reciprocally as a straight and rigid rod, and therefore generates no thrust (scallop theorem). In the limit Sp 1, viscous forces dominate and the amplitude of the bending wave declines rapidly as it travels along the length of the filament. The rest of the filament remains essentially straight, not generating thrust and impeding locomotion due to drag. Analytical and numerical solutions therefore predict the existence of an optimum propeller for an intermediate value of Sp [16,19,20,17,21]. This prediction has been confirmed experimentally [22,23]. Yu et al. showed that data for different filament lengths and diameters all collapse onto a single curve of propeller performance (measured in terms of normalized thrust force) vs. Sp for a given propeller design [23]. A crucial implication of this is that the relationship between propeller performance and Sp is fixed with respect to the parameters that Sp accounts for -namely the geometry (L and r), actuation frequency (ω), and material properties of filament (E) and fluid (µ). In other words, as long as one has the ability to control at least one of these parameters, it is possible to tune Sp to achieve optimum performance with respect to these parameters [17,23].
In our work, E, µ, and ω were given by the elastic modulus of PDMS, viscosity of Percollculture media mixture, and the contraction frequency of the neuromuscular unit, respectively. Cross-sectional dimensions (thereby the equivalent r) were bounded within a narrow range primarily by fabrication constraints, leaving the filament length L as the most feasible tuning parameter. For the double-tail swimmer design, our simulations for different tail thicknesses and lengths predicted a maximum swimming velocity of approximately 0.96 µm/s. For each tail thickness, this maximum velocity was reached at a different tail length (Fig. S4). Experimental results obtained for 10 µm tail thickness were in close agreement with the simulation results ( Fig. 4E in main text). Remarkably, when we plot our model predictions and experimental results against Sp instead of tail length, the results show that the three different curves for different tail thicknesses converge to the same trend. Maximum velocity (0.96 µm/s in simulation, 0.85 µm/s in experiment) is achieved for Sp ≈ 2. Velocity declines to zero for smaller Sp and plateaus at a suboptimal value for larger Sp (Fig. S10), consistent with theoretical predictions [16,19,20,17,21]. We therefore conclude that the performance of our double-tail swimmer design is optimized with respect to tail geometry, material properties, and actuation frequency -for the given actuation design (two tails actuated cyclically at their base with a displacement amplitude of 20 µm).
For further optimization of thrust, we must consider factors that the Sperm number analysis does not take into account. One such factor that we have examined in this study is the number of propellers and their arrangement. Our models predicted that adding a second tail improved swimming speed by an order of magnitude ( Fig. 3D and Fig. S4). While the additional thrust force from the second tail undoubtedly contributes to improved performance, we believe that the major advantage of the double-tail design is the symmetrical arrangement of the tails which ensures that the swimmer's head doesn't swing laterally, as it does in the single-tail design (Movie S5). This lateral swinging of the head is a major energy drain due to the large viscous drag on the head, and may therefore be impeding swimming performance. Improved performance of a double-tailed swimmer over a single-tail design was also theoretically predicted by Lauga, who likewise postulated that the reason was the elimination of lateral swinging of the swimmer head [17]. It is conceivable that adding more tails would improve swimming velocity, as long as an even number of tails are used and are arranged symmetrically.
Another viable strategy for further improvement of propulsion force is to increase actuation amplitude (Sp accounts for the actuation frequency but not the amplitude). In our experiments, the maximum active contraction force exerted by the muscle strips on the scaffold legs was approximately 7.64 µN (F peak in Eq. S6), corresponding to 20 µm displacement at the base of each tail (leg deflection in Fig. 3B and Fig. 4D). While this force is already significantly higher than the values reported by previous studies using the same cell source, similar tissue size, and scaffold rigidity (∼ 2 µN in [24] and ∼ 1.5 µN in [1]), there are several possible strategies to improve it further. Skeletal muscle differentiation efficiency and force generation can be enhanced by treating the cells with IGF-1 (insulin-like growth factor-1) [25,26]. The force output to the scaffold legs (F out in Eq. S5) can also be improved by refining tissue culture techniques to increase muscle fiber alignment and the ratio of the volume occupied by contractile fibers to the volume of passive material (i.e. non-contractile cells and ECM) [27]. The displacement generated due to muscle contraction also depends on the stiffness (i.e. the spring constant, k) of the scaffold legs. However, changing scaffold stiffness will also affect muscle forces. The displacement-stiffness relationship for a muscle-based actuator is therefore nonlinear. Myogenic differentiation and formation of contractile fibers requires tension in the tissue [28], which cannot be developed if the scaffold is too soft. We therefore expect muscle force and resulting leg deflection to diminish in the limit k → 0. On the other hand, muscles cannot generate arbitrarily large forces, so the leg deflections will also diminish in the limit k → ∞, indicating the existence of an optimum leg stiffness. Maximizing actuation amplitude therefore requires careful characterization of the displacement-stiffness relationship for a muscle-based actuator. To estimate how much swimming velocity may be improved by increasing muscle actuation output, we performed simulations for "hypothetical actuators" by artificially changing the value of muscle generated force in simulation to obtain different leg deflection values within a reasonably realistic range of 10 and 40 µm (our current actuator gives 20 µm). The results illustrate that it may be possible to attain several times higher swimming velocities by increasing actuation amplitude within this range (Fig. S11).
We note that in our swimmer design, the thin PDMS tails are the predominant source of thrust, and the contribution of the rest of the body to thrust force is negligible: While the legs of our swimmer scaffold are flexible enough to deform under muscle forces, they are too rigid for bending waves to propagate along them. To give a comparison, the thickness (h) of the tails and legs of our swimmer are 10 µm and 100 µm, respectively. Since the moment of inertia (I) scales as h 3 , bending force for the legs (∼ EI/L 3 per unit length) is 3 orders of magnitude higher than that for the tails. Hence, the legs operate at the Sp 1 limit and do not generate significant thrust due to the scallop theorem. The head of the scaffold is even more rigid and does not deform appreciably. We therefore evaluate the influence of the size of the body (head-legs-tissue assembly) on swimming performance only in terms of its contribution to drag. In the low Re regime, the drag force on an object is proportional to the body size (see Eq. S9 and corresponding discussion). Hence, reducing the size of the swimmer body for the given tail actuation will improve swimming velocity. Indeed, theoretical analysis has predicted significantly higher swimming velocities for swimmers with slender heads compared to swimmers with spherical heads [17]. For our swimmer design, the limiting factor for miniaturization is the neuromuscular tissue which determines the size and shape of the scaffold. Thinner and shorter muscle strips could be engineered by refining tissue formation and culture procedures [24,29], and neurons could conceivably be cultured as an adherent monolayer on the scaffold surface instead of being seeded as a large spherical cluster. It is important to note, however, that these modifications will most likely affect the muscle contraction amplitude as well as the drag. Minimizing drag by miniaturization of the swimmer body therefore requires careful investigation of the relationship between contraction force and the size of the muscle tissue.
Lastly, we consider the implications of increasing swimming velocity. For a given body length, higher velocity will correspond to a higher Reynolds number. Since the theoretical analyses on which we based our discussion of swimmer performance are valid strictly for Re 1, the optimization strategies we have provided here may no longer apply for Re ∼ 1 and higher. Due to the fundamental differences in the physical character of low Re and high Re flows, swimming at these two regimes typically involves distinctly different propulsion mechanisms. While flagellar and ciliary propulsion is common at low Re, swimmers at high Re employ mechanisms such as the flapping or undulation of the entire body or large fins [30]. Therefore, when building a biohybrid swimmer, the targeted flow regime and an appropriate propulsion mechanism should be determined at the beginning of the design process.       shows that performance is recovered, but sustained for only about 8s, followed by fatigue and deceleration. (Right column) After allowing a longer 10min rest period, steady-state swimming performance is recovered and sustained (albeit intermittently) for approximately 30s before beginning to decelerate. In total, the swimmer travels a similar distance (25µm) during the third recording as it did during the first.     [32] 10 kPa Muscle density [33] 1.06 g/cm 3 Dynamic viscosity [34] 1.15 mPa.s

Movie Legends
Movie S1. Tail deflections induced by neuromuscular actuation. Before releasing the swimmer, tails are cut first and the capability of muscle contractions to induce tail deflection is verified. Movie plays at 2x speed.
Movie S2. Vertical distance between swimmer and bottom surface. Z-stack taken after suspending the swimmer in Percoll-media mixture, covering a range starting below the bottom surface of the dish up to a plane above the swimmer. Movie S5. Simulation of the single-tail swimmer design illustrating the lateral swinging motion of the swimmer head.