Fluid Dynamics of Ballistic Strategies in Nematocyst Firing

Abstract

Nematocysts are stinging organelles used by members of the phylum Cnidaria (e.g., jellyfish, anemones, hydrozoans) for a variety of important functions including capturing prey and defense. Nematocysts are the fastest-known accelerating structures in the animal world. The small scale (microns) coupled with rapid acceleration (in excess of 5 million g) present significant challenges in imaging that prevent detailed descriptions of their kinematics. The immersed boundary method was used to numerically simulate the dynamics of a barb-like structure accelerating a short distance across Reynolds numbers ranging from 0.9–900 towards a passive elastic target in two dimensions. Results indicate that acceleration followed by coasting at lower Reynolds numbers is not sufficient for a nematocyst to reach its target. The nematocyst’s barb-like projectile requires high accelerations in order to transition to the inertial regime and overcome the viscous damping effects normally encountered at small cellular scales. The longer the barb is in the inertial regime, the higher the final velocity of the projectile when it touches its target. We find the size of the target prey does not dramatically affect the barb’s approach for large enough values of the Reynolds number, however longer barbs are able to accelerate a larger amount of surrounding fluid, which in turn allows the barb to remain in the inertial regime for a longer period of time. Since the final velocity is proportional to the force available for piercing the membrane of the prey, high accelerations that allow the system to persist in the inertial regime have implications for the nematocyst’s ability to puncture surfaces such as cellular membranes or even crustacean cuticle.

1. Introduction

Nematocysts are specialized organelles possessed by members of the phylum Cnidaria (e.g., jellyfish, corals, anemones) and certain microorganisms such as dinoflagellates [1,2]. Nematocysts are used by organisms for prey capture and feeding, but also for defense, transport, digestion and other various functions [3,4]. There is a broad diversity in the morphology and function of nematocysts, with about 30 different types classified, and multiple types often found in the same organism [5,6]. However, each type has a similar basic structure comprising a fluid-filled capsule containing a highly-coiled tubule tipped with a barb-like stylus. Schematics of several nematocysts are shown in Figure 1. Our study focuses on a representative type of nematocyst commonly found in Hydra and the subphylum Medusozoa: the stenotele [7] (illustrated in Figure 1D and Figure 2). Similar nematocysts are found in dinoflagellates [8]. When triggered, the capsule ejects the stylus which pierces its target, allowing the venom-carrying tubule to be injected. The ejection process for a stenotele nematocyst is shown in Figure 2. Nematocysts are of interest as a model for applications such as microinjector design for targeted drug delivery [9,10], and there has been much interest in categorizing and carefully classifying the characteristics and structures of each component of the nematocyst in order to better describe and quantify the firing mechanism, which consists of a lid (operculum) opening followed by tubule discharge [6,11].

While the speeds involved are modest across metazoa, averaging in the range of 9–18 m/s [12], the ejection of the stylus during nematocyst firing is claimed to be one of the fastest-known accelerations in the animal world at any scale, generating accelerations in excess of 5,000,000 g [12,13]. Moreover, the speeds are relatively high for a microscale process. The upper limit of the acceleration is not precisely known because the time scales involved (microseconds) are beyond what can currently be resolved using high speed videography [12,14]. In some cases, mathematical models have been used to analyze the kinematics of the nematocyst firing [15], but the role of fluid dynamic interactions and boundary layer effects in small-scale ballistics are not well understood [8,16]. Previous mathematical modeling of general biological ballistics involved computing the velocity of the projectile from a force balance equation comprising a motor, spring, and latch without taking the fluid dynamics of the medium into account [17] or by using simplified fluid approximations [18].

Cnidarians are obligately aquatic, and we hypothesize the interaction of water with the stylus of a nematocyst significantly affects the motion of the projectile as well as the ability to reach the prey. The physics governing the motion of structures in water at the microscale are usually over-damped because the fluid environment is dominated by viscous forces rather than inertial forces [19]. The Reynolds number, $Re$, is a nondimensional number that quantifies the ratio of inertial to viscous forces in a fluid. When the characteristic length scale of a system is relatively small, such as for unicellular organisms, 1.0–10.0 ×${10}^{-5}$ m, then $Re\ll 1$, and viscous forces dominate the fluid flow. As such, most cellular processes and cellular motility occur in highly viscous regimes. The fluid dynamics of such cellular environments can be non-intuitive.

For example, organisms moving in these micro-scale environments at typical speeds (0–0.007 m/s, [20]) cannot achieve net displacement if reciprocal motions are used for propulsion. This result is known as the the scallop theorem [19]. Only non-reciprocal motion, such as the rotation of a flagellum of a bacterium, can result in net movement. Furthermore, boundary layer effects are significant. The effective boundary layer of a cell may be more than ten times its radius. Note, however, that nematocyst firing is in the intermediate $Re$ range at $Re\sim 900$ where inertial forces are also significant. We expect that the fluid dynamics of this system is fundamentally different from most microscale systems.

In our previous study, [21], we examined the fluid dynamics resulting from firing a micro-sized nematocyst stylus/tubule complex into a viscous fluid. From estimates inferred from experimental measurements, the acceleration of the barb results in a transition from the viscous regime to the inertial regime ($Re\ll 1$ to $Re\sim 900$). The simulations comprised a ballistic structure (hereafter called the barb) and prey, with both structures immersed in viscous fluid governed by the incompressible Navier-Stokes equations. The barb was accelerated to a final constant velocity using tethered forcing points, where discrete points compromising the barb are physically linked via stiff elastic springs to a corresponding collection of points that move with a specified velocity (described in Section 2.4). Note that after the acceleration phase, a constant preferred velocity of 18 m/s was maintained using these target points. In this paper, the barb is accelerated to a final velocity as before. However, after the short acceleration period, the attachment to the tether forces is removed, and the barb’s motion is due only to the interaction with the fluid. Said another way, the barb coasts to reach the prey.

Results from [21] showed that the barb had to overcome the viscous regime in order to reach the prey due to the size of the boundary layer that had the effect of pushing the prey out of the way. At small $Re$, these large boundary layers around the barb resulted in a front extending ahead of the barb. Our results showed that this effect where the prey was pushed out of the way was present regardless of barb length. At higher $Re$, the barb was able to reach the target due to the reduced boundary layer size. Simulation data also showed the barb more easily reached the larger prey at biologically relevant parameter values reported for Hydra spp. We found the momentum imparted by the larger additional mass from the surrounding and enclosed fluid of the prey, which resulted in increased resistance to being pushed out of the way by the moving fluid front generated by the barb. In contrast, results from model simulations in this work show that the size of the prey does not affect the barb’s ability to reach the prey at the biologically relevant value of $Re=900$. We find the size of the barb contributes significantly to the average and final distance between the prey and the barb (data presented in Section 3.4 and discussed in Section 4).

The use of numerical simulations of nematocyst discharge is a powerful tool for understanding the fluid dynamics involved given that the fast kinematics prohibit direct experimental measurement. Our previous work did not, however, consider the deceleration of the barb after the initial acceleration period. In our current study we focus on the barb’s initial approach following the launch towards the prey where the barb is allowed to coast through the fluid toward its target. This approach allows us to determine the effects of the osmotic pressure gradient that drives the barb during firing [15] and the subsequent interaction of the boundary layers of the projectile and prey that act to both slow the barb and push the prey out of the way. By varying the $Re$, prey size, and barb length, we determine the necessary conditions on the initial acceleration of the ballistic structure in order to reach the prey with a non-negligible final speed. Similar to our previous study, two-dimensional immersed boundary simulations are used to allow the exploration of a wide parameter space to identify key phenomenon. As these simulations require us to solve the fully-coupled fluid-structure interaction problem using the full Navier-Stokes equations at high resolution, it is not currently feasible to perform a large number of three-dimensional simulations.

Of particular interest in the current study is how the transition to the inertial regime and the associated changes in boundary layers affect the speed of approach and ultimately whether or not the barb contacts the prey. Importantly, the large magnitude accelerations involved in nematocyst firing result in a changing fluid environment and a more complicated interaction than those seen in strictly lower $Re$ microscale systems [21,22,23,24]. Because the types and functions of nematocysts are highly varied, our study focuses on the general structure of the basic interactions rather than trying to resolve the fine structures and mechanisms of a particular class of organelle. In this way, we seek to understand the underlying fluid dynamics of interactions across changing Reynolds numbers.

2. Materials and Methods

2.1. Model Overview

The idealized prey capture model consists of a relatively stiff barb and a flexible, circular prey immersed in a viscous fluid (water). Both the barb and the prey resist bending and stretching. A schematic of model components is shown in Figure 3. Each component of the model is described in detail below. The domain was taken to be a square with length and width 10 times the base barb length. The barb and prey were centered horizontally in the domain and evenly spaced vertically with the midpoint between them at the origin.

For the dimensions of the barb and prey, we chose parameters based on typical biological measurements [7,8,11,12] (

Table 1. Parameters taken from measurements for organisms. The prey radius can be much larger, for example in Hydra. Experimental values for the barb size in Hydra are difficult to obtain because the tube shaft is still attached to the barb during the ejection process and the ultrafast dynamics of firing. Note the length of the barb and tubule together can vary over several orders of magnitude, for details on the diversity of structures and prey sizes, see the references noted in the table. We consider a typical structure size on the order of the distance between the structures to focus on the contact event.

Typical Sizes for Nematocysts and Prey
Structure	Length	Source
Length of nematocyst tubule and barb	3.0– > 12.0 × ${10}^{-5}$ m	[12,25]
Spacing between objects	1.0–5.0 × ${10}^{-5}$ m	[7,8,11,12]
Radius of prey	0.5–2.5 × ${10}^{-5}$ m	[12]

). The parameters used in the computational model are summarized in

Table 2. Parameters used in simulations. The size of the fluid domain was chosen to be 10 times the length of the structures. The spacing between the length of the barb and the prey was chosen to be the same length as the barb in the base case so that the barb had to travel at least one body length to reach the target prey. * The Reynolds number is calculated using the maximum achieved barb velocity at the end of the acceleration phase. * The mesh size varies from 128–512 points in both the x and y directions, and $dx$ is reported for the finest mesh of the adaptive mesh grid.

Simulation Parameter Values
Structure	Length
Length of barb (L)	$5.0\times {10}^{-5}$ m
Spacing between objects (h)	$5.0\times {10}^{-5}$ m
Radius of prey (R)	$2.5\times {10}^{-5}$ m
$Re$	0.9–900 *
Domain size ($\mathsf{\Omega }$)	$5.0\times {10}^{-4}$ m
$dx$ *	$9.766\times {10}^{-7}$ m
$k_{spring}$	$1.64\times {10}^{11}$ N/m
$k_{beam}$	$7.65\times {10}^{10}$ N·m
$k_{target}$	$1.64\times {10}^{11}$ N/m

. Note that the initial distance between the barb and the prey is equal to the length of the barb in the base case so that the barb must travel at least one body length to reach the prey. The rationale for this choice of distance is based on the fact that the barb which punctures the prey is attached to the organelle. It is reasonable to assume that the barb will not travel great distances in order to reach the prey, and in many instances, such as for a stenotele nematocyst, the barb everts as part of the firing mechanism so traversing a distance of a body length is a reasonable starting value.

2.2. Barb Construction

The base length of the barb is set to $L=5.0\times {10}^{-5}$ m, a typical size for Hydra. Boundary points representing the barb are directly linked by springs to a set of tether points (described in more detail below) moving at a specified velocity. The tether and barb points are initialized at the same location. As the preferred position of the barb changes in time, a target force is applied to the boundary which accelerates the barb forward. This force is removed during the coasting phase (see details below).

2.3. Prey Construction

The radius of the circular prey is set to $R=2.5\times {10}^{-5}$ m, a typical prey size for Hydra (Table 2). The structure representing the prey is neutrally buoyant and untethered, representing a suspended single-celled organism that is free to deform. Our primary interest is the fluid dynamics of the interacting boundary layers, as such, we do not consider any complex geometry or material properties of the prey. The circular model simply resists bending and stretching. The stiffnesses were chosen such that small deformations are observed during the simulation.

2.4. Movement of the Barb

In order to accelerate the barb from rest to its maximum velocity, tether points that impose a preferred trajectory are used to drive the barb during the acceleration phase. Tether forces are a common way of imposing active motion in IB simulations [26,27]. During the acceleration period, the force density along the length of the barb is calculated by

$$\begin{array}{c}\hfill {\mathbf{F}}_{target}=k_{target}({\mathbf{X}}_{barb}(s,t)-\mathbf{Y}(s,t))\end{array}$$

(1)

where $\mathbf{Y}(s,t)$ is the preferred (target) position of the material point whose current position is ${\mathbf{X}}_{barb}(s,t)$. Both the target point and actual boundary points of the barb are initialized at the same locations. As the simulation evolves, differences between the preferred and actual positions of the material points impose “penalty” forces at the material point locations, driving the material points towards the preferred positions. The stiffness of these tether springs, $k_{target}$, is chosen so that the distance between the target point and its corresponding material point is not more than $1\%$ of the barb’s length throughout the simulation (see Table 2 for exact value of $k_{target}$).

During the initial acceleration phase, the barb is accelerated over a time period of $2.78\times {10}^{-7}$ s to a final velocity of 18 m/s, resulting in a constant acceleration of $6.4\times {10}^7$ m/s${}^2$ from $t=0$ to $t=2.78\times {10}^{-7}$ s. The maximum velocity and acceleration times are in the ranges estimated from [8,15]. After this initial acceleration, the barb is allowed to freely coast without influence from the tether points for an additional $7.72\times {10}^{-6}$ s for a total simulation time of $8\times {10}^{-6}$ s, representing a biologically relevant estimate for both the time of acceleration and the time for the barb to contact the prey [8,15]. This coasting phase is achieved by setting the stiffness constants of the springs connecting the barb to the tether points are set to zero.

Table 3. Typical values for parameters relating to acceleration of the barb and time stepping.

Acceleration and Time Parameter Values
Time duration	$8\times {10}^{-6}$ s
Maximum barb velocity	18 m/s
Ramp up time	$2.78\times {10}^{-7}$ s
Total number of time steps	2,770,083
Number of time steps during initialization	96,153
Simulation time	$8\times {10}^{-6}$ s
Time step ($\mathsf{\Delta }t$)	$2.888\times {10}^{-12}$ s

shows the temporal parameters used in the simulations.

2.5. Material Properties of Barb and Prey

For both the prey and the barb, the structures are modeled to resist stretching and bending so that the total force on each is given by the equations

$$\begin{array}{cc}\hfill {\mathbf{F}}_{barb}& ={\mathbf{F}}_{target}(s,t)+{\mathbf{F}}_{spring}^{barb}(s,t)+{\mathbf{F}}_{beam}^{barb}(s,t),\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill & \hfill {\mathbf{F}}_{prey}& ={\mathbf{F}}_{spring}^{prey}(s,t)+{\mathbf{F}}_{beam}^{prey}(s,t),\hfill \end{array}$$

(3)

where ${\mathbf{F}}_{spring}(s,t)$ represents the restoring force from the springs resisting stretching and ${\mathbf{F}}_{beam}^{prey}(s,t)$ are forces due to bending resistance. In the immersed boundary (IB) framework, tension in structures can be modeled by

$$\begin{array}{c}\hfill T(s,t)=k_{spring}\left(\left|\frac{\partial {\mathbf{X}}_{prey}}{\partial s}\right|-1\right)\end{array}$$

(4)

where $k_{spring}$ is the spring stiffness associated with the given structure. The forces due to stretching can then be calculated by

$$\begin{array}{c}\hfill \mathbf{F}{(s,t)}_{spring}=\frac{\partial }{\partial s}\left(T\tau \right),\end{array}$$

(5)

where $\tau$ represents the unit tangent vector of the parameterized curve $\mathsf{\Gamma }$ describing the structure

$$\begin{array}{c}\hfill \tau =\frac{\partial {\mathbf{X}}_{prey}/\partial s}{\left|\partial {\mathbf{X}}_{prey}/\partial s\right|}.\end{array}$$

(6)

To model the bending forces, we follow the work of [28] and discretizing the energy functional

$$\begin{array}{c}\hfill E_{bend}=\frac{k_{bend}}{2}{\int }_{\mathsf{\Gamma }}{\left|\frac{{\partial }^2\mathbf{X}(s,t)}{\partial s^2}\right|}^{2}ds\end{array}$$

(7)

where $k_{bend}$ is the bending rigidity. Note that this force is calculated using the Fréchet derivative [29]. The preferred curvature for both the barb and the prey is taken to be zero along the entire boundary. Values for the parameters $k_{spring}$ and $k_{bend}$ are located in Table 2.

2.6. Immersed Boundary Method

The fully-coupled fluid-structure interaction model describing the moving barb, the passive, flexible prey and the surrounding fluid is numerically solved by implementing the model in an immersed boundary framework [29,30]. The immersed boundary method has been successfully used to model a variety of biological fluid dynamics problems at intermediate $Re$ ($0.01<Re<1000)$ (e.g., [31,32,33]). The elastic structures describing the barb and the prey are defined in a moving Lagrangian framework, while the surrounding fluid is described in a fixed Eulerian framework using the Navier-Stokes equations for a viscous, incompressible fluid,

$$\begin{array}{ccc}\hfill \rho ({\mathbf{u}}_t(\mathbf{x},t)+\mathbf{u}(\mathbf{x},t)·\mathbf{\nabla }\mathbf{u}(\mathbf{x},t))& =& -\mathbf{\nabla }p(\mathbf{x},t)+\mu {\mathbf{\nabla }}^2\mathbf{u}(\mathbf{x},t)+\mathbf{f}(\mathbf{x},t),\hfill \end{array}$$

(8)

$$\begin{array}{cccc}\hfill & \hfill \mathbf{\nabla }·\mathbf{u}(\mathbf{x},t)& =& 0,\hfill \end{array}$$

(9)

where $\mathbf{u}$ represents the fluid velocity $\mu$ represents the dynamic viscosity of the fluid, p represents the fluid pressure, $\rho$ represents the constant fluid density, and $\mathbf{f}$ represents the forces contributed by the immersed structures. Dividing Equation (8) by $\rho$ and non-dimensionalizing yields

$$\begin{array}{ccc}\hfill Re\left(\frac{\partial \mathbf{u}}{\partial t}+\left(\mathbf{u}·\mathbf{\nabla }\right)\mathbf{u}\right)& =& \mathsf{\Delta }\mathbf{u}-\mathbf{\nabla }p+\mathbf{f},\hfill \end{array}$$

(10)

where $Re$ is the Reynolds number with characteristic velocity U taken to be the velocity of the barb at the end of the acceleration phase and with the characteristic length L taken to be the length of the barb. $Re$ is altered in the system by changing $\mu$, which allows us to only vary one parameter in the Navier-Stokes equations while leaving the length, elastic properties of the structures, and time scales unchanged. The external forcing term $\mathbf{f}$ is calculated by spreading the singular force density located along the Lagrangian structure $\mathbf{F}$ according to the equation

$$\begin{array}{c}\hfill \mathbf{f}\left(\mathbf{x},t\right)={\int }_{\mathsf{\Gamma }}\mathbf{F}(s,t)\delta (\mathbf{x}-\mathbf{X}(s,t))ds,\end{array}$$

(11)

where $s\in \mathsf{\Gamma }$ are the material points, $\mathbf{X}(s,t)$ is the position of the material point s at time t, and $\delta (\mathbf{x},t)$ is the Dirac delta function. To determine the velocity of the Lagrangian points in time, the local fluid velocity is interpolated using the equation

$$\begin{array}{c}\hfill \mathbf{U}\left(s,t\right)={\int }_{\mathsf{\Omega }}\mathbf{u}(\mathbf{x},t)\delta (\mathbf{x}-\mathbf{X}(s,t))d\mathbf{x},\end{array}$$

(12)

where $\mathsf{\Omega }$ is the fluid domain, ensuring the system observes the no-slip condition for a viscous incompressible fluid.

2.7. Numerical Methods

The base case for the barb structure was discretized using 100 evenly spaced points, corresponding to a spacing of $ds=dx/2$, where $dx$ is the spacing of the underlying fluid grid at the finest level of an adaptive grid. The number of points for the other barb cases was scaled proportional to the length to preserve the $ds$ spacing. The prey structure was initialized as a closed circle with points evenly spaced $ds$ apart as

$$\begin{array}{c}\hfill {\mathbf{X}}_{prey}=R\left(cos(s/R),sin(s/R)\right),0\le s\le 2\pi R.\end{array}$$

(13)

To numerically solve the IB equations, we implemented our model within the open-source software distribution IBAMR version 0.5.0, which uses a formally second-order accurate solver in time and space, with decreased accuracy for discontinuous functions including pressure jumps across thin interfaces. IBAMR is a parallelized distribution with adaptive-mesh refinement, which allows high refinement near the boundary structures and in regions of high vorticity with relatively low computational cost in the far-field fluid calculations [34]. For our simulations, the grid spacing varies dynamically across up to four levels of refinement so that $dx$ varies from $\mathsf{\Omega }/128-\mathsf{\Omega }/512$, where $\mathsf{\Omega }$ is the domain size. Using these values, finest grid spacing is $dx=9.77\times {10}^{-7}$ m (see Table 2). Implementing the current models in an IB framework with adaptive mesh refinement facilitates the building of more computationally challenging models in future work, such as multiple nematocysts firing in sequence as well as extensions of the model into three dimensions.

3. Results

We first present the results of the biologically relevant base case, with $Re\sim 900$, the length of the barb $L=5.0\times {10}^{-5}$ m and the radius of the prey $R={10}^{-5}$ m. Other parameter values are listed in Table 2. Recall that the barb is initially driven at a constant acceleration and is then released to coast towards the prey. A time series of the vorticity and velocity vectors is shown in Figure 4. The boundary layer around the barb is visible by inspecting the color map of the vorticity. As the barb approaches the prey, the shearing of the fluid introduces a layer of non-negligible vorticity around the barb. This boundary layer may also be observed through the deflections in the local velocity vectors. The prey is slightly pushed upward, but minimal deformations of the boundary or shearing of the fluid around the prey is observed. The barb has contacted the prey by the fourth frame (bottom right) of Figure 4. These results are similar to those found for a barb moving at a constant prescribed velocity at $Re=900$ reported in [21].

Inspection of this time series shows that the vorticity around the barb is strongest during the acceleration period. As the simulation progresses, the speed of the barb decreases as it approaches the prey. Note that this is significant since the speed of the barb at contact will be proportional to the puncture force (more details below). The interaction of the boundary layers around the barb and the prey serves to both slow the barb down and to push the prey slightly upwards. At the biologically relevant $Re$ of 900, the barb has enough inertia to easily reach the prey.

3.1. Reynolds Number

For the results described below, $Re$ is varied by changing the dynamic viscosity $\mu$. Viscosity alone was varied so that all length and time scales could be fixed. Given that force does not scale with velocity at higher $Re$, it was also possible to fix the elastic properties of the prey and barb. Figure 5 shows the vorticity and velocity of the fluid at $Re=0.9,9,900$ after 3/4 of the total simulation time has elapsed. Because the barb is coasting with significant deceleration, the effective $Re$ is changing throughout the simulation. For ease of comparison, the $Re$ is calculated using the maximum velocity of the barb just before coasting.

Significant changes in the boundary layers can be seen across these $Re$. In general, as the Reynolds number decreases, the thickness of the boundary layer increases and the slow-down effect is most pronounced. In the most viscous case $Re=0.9$, there is significant fluid velocity and vorticity far from the prey and barb. Given the large volume of fluid entrained by the barb, the prey is pushed forward and out of the way, and the barb rapidly decelerates. Note that the prey in the snapshot has been pushed slightly upwards and this displacement continues during the rest of the simulation. At $Re=9$, the volume of fluid entrained is significantly reduced, through there is still non-negligible flow far from the barb. This region of entrained fluid is then further reduced at the higher value of $Re=900$. The prey is displaced upwards once contact occurs (not shown).

Figure 6 and Figure 7 show temporal snapshots of the color maps of the vorticity and the velocity vectors for simulations performed at $Re=0.9$ and 9, respectively. For $Re=0.9$, the barb rapidly entrains a large volume of fluid, and oppositely spinning vortices are observed on either side of the barb. Note that this flow pattern is not observed at the higher $Re$. As a consequence of the large boundary layer and the formation of these vortices, the prey is pushed upward throughout the simulation. For $Re=9$, vorticity around the barb is more localized and stronger. The size of the boundary layer is non-negligible, and the barb begins to push the prey upwards during its approach.

To better understand the effects of $Re$, we compare the distances traveled by each barb after 3/4 of the simulation has completed. In previous work, the motion of the barb was controlled by tether points driving it at the maximum achieved velocity throughout the simulation. In this case, the barb travels about 1.5 body lengths after 3/4 of the total simulation time had elapsed. In all of the current simulations, the barb has not reached the prey after 3/4 of the time has elapsed but has traveled $0.95L$, where L represents the length of the barb (see Table 2). When $Re=9$, the barb travels about $0.72L$ after 3/4 of the simulation, and at $Re$ = 0.9, the barb travels $0.55L$. These significant slow downs are due to the increasing effects of the boundary layer interfering with the forward momentum of the barb.

The distance between the barb and prey over time is shown in Figure 8 for $Re=0.9,9$, and 900. We set a threshold of 4 grid spaces to indicate that the barb touches the prey. Since the grid spacing $dx$ varies during a simulation, we use the most refined grid level to compute the touch distance, $dx=\mathsf{\Omega }/512$ (see Table 2). This choice is due to the fact that regularized immersed boundary points have an effective radius of about 1.5 $dx$, and we must account for the thickness of both the barb and the prey [35]. We note that at this distance both the barb and prey appear to be ‘connected’ and travel together.

In Figure 8, we observe that at any point in time during the simulation, the distance between the barb and prey monotonically decreases as $Re$ increases. Furthermore, the distance between the barb and the prey decreases monotonically as the simulation time increases for all $Re$. For the case of $Re=0.9$, the distance decreases only about 20%, and the graph is nearly linear. When $Re=9$, the distance between the barb and prey decreases about 75% during the simulation, and a nonlinear relationship between distance and time is observed. For $Re=900$, the distance initially rapidly decreases and then plateaus. The flat portion observed at the end of the simulation is likely due to the effective width of the regularized boundary layer.

3.2. Barb Length

We vary the barb length at $Re=900$, keeping the values for bending rigidity and elastic stiffness of the structures the same. Figure 9 shows temporal snapshots of a color map of the vorticity and the velocity vector fields during a simulation for a barb whose length is 2.5 times the base value. The region of vorticity around the barb is narrower and of stronger vorticity (compare to Figure 4). The barb and prey are effectively in contact and move together at $t=5.5\times {10}^{-6}$ s.

The following trends are observed from simulations where the barb length is varied from $L=5.0\times {10}^{-5}$ m. For barbs with lengths equal to or greater than the base case, the barb experiences a noticeable slow-down as it moves closer to the prey, but it is still able to reach the prey. The longer barbs reach the prey more quickly than when $L=5.0\times {10}^{-5}$ m. Barbs of smaller lengths exhibit wider boundary layers (data not shown). As a consequence, the relatively larger volumes of entrained fluid slow the short barbs down long before they reach the prey. These effects are more pronounced at lower $Re$. For lowest $Re=0.9$ simulations, the longest barbs are still unable to contact the prey. Data for simulations that reach the prey over a range of different barb lengths and $Re$ are presented in Section 3.4.

Figure 10 shows the distance between the barb and prey over time during a simulation for several different values of barb length when $Re=900$. At any given time during the simulation, the distances between the barb and prey decrease with increasing barb length. This effect may be explained by noting that the region of fluid sheared is narrower for longer barbs, indicating thinner boundary layers. For these longer barbs, the entrained volume of fluid per unit length of the barb is consequently lower, and longer barbs are able to coast for longer periods of time. Furthermore, the boundary layer is smaller such that the prey are not pushed out of the way.

3.3. Prey Size

For each $Re$ we varied the prey size from a radius of 1/4 the base case to 4 times the base case. The prey size does not appreciably change the distance between the prey and the barb for simulations at $Re=900$, the biologically relevant case. The only noticeable difference at $Re=900$ is that the larger prey would resist being pushed by the barb’s approach, and the smallest prey would move upwards, likely due to the smaller overall contained volume of fluid in the closed prey structure. Prey size did have an effect on the distance between the barb and prey for $Re=9,0.9$, and this is discussed in more detail below.

3.4. Barb Length Determines Whether the Barb Reaches the Prey for Intermediate Values of $Re$

The average distance between the prey and barb over the course of a simulation as well as the distance at the end of the simulation was computed for several values of barb length and prey radius (data shown in Figure 11). Distance values less than or equal to the threshold value of 4 grid spaces were set to zero for the data in Figure 11. The justification for our choice of the threshold value is provided in Section 3.1. The average distance was computed during the entire coasting period, from the time the barb was released to the end of the simulation, i.e., for time values from $6.48\times {10}^{-7}\le t\le 8\times {10}^{-6}$. The graphs of the final distance (Figure 11B) display the distance between the barb and the prey at the end of a simulation. Note that the only parameter values when barb hit the prey by the end of a simulation occur when $Re=900$ and $L\ge 5.0\times {10}^{-5}$ m (bottom right panel of Figure 11). The average distance is indicative of how quickly the barb approaches the prey during the course of a simulation.

We observe a general trend that as $Re$ increases, the average distance between the prey and barb decreases. The barb reaches the prey only for data corresponding to zero values for the final distance when $Re=900$. For the case of $Re=0.9$, the final distance is greater than $1.1\times {10}^{-5}$ m from the prey for a range of values for barb length and prey radius. For the smallest value of $Re$, the average distance decreases as the prey radii increases. As the $Re$ increases, the average distance becomes less sensitive to the size of the prey. For $Re=900$, the average distance between the barb and the prey depends on the size of the barb and is insensitive to the radius of the prey, as illustrated by the horizontal contour lines for the distance data when $Re=900$ in the bottom row of Figure 11.

4. Discussion and Conclusions

The results from the simulations show that in order for the barb to reach the prey with velocities sufficient to potentially puncture the prey, large velocities due to high initial accelerations are needed to overcome the strong effects of the boundary layer interactions. At lower $Re$, large boundary layers result in dissipation of the momentum gained during the rapid acceleration phase. At the lowest $Re$ examined, the barb stalls after traveling less than half its body length. Our results also suggest that the acceleration of the added mass of the surrounding fluid is important. For longer barbs, narrower boundary layers develop leading to a smaller volume of entrained fluid per unit length of barb. As such, the length of the barb can play an important role in the barb’s ability to reach and penetrate its target. Overall, this study indicates that the boundary layer interactions involved in decelerating the barb and potentially pushing the prey out of the way are quite complicated due to the rapid transitions through fluid regimes in several parts of the firing event.

These results are a departure from our previous study where the barb moved with a constant velocity for most of the simulation. Data from simulations using the earlier model indicated that the size of the prey, rather than the size of the barb, was more influential in the interactions [21]. This is likely due to the added momentum imparted from the target points to the barb, sustaining the motion of the barb and overcoming the boundary layer effects. In the previous work, the larger concern was the buffeting of the prey away from the barb as the boundary layers interacted, particularly at lower $Re$. In our current simulations, the same buffeting effects at lower $Re$ occur but were secondary to the slow-down effects of the barb, indicating that the barb requires sustained external forcing to maintain its velocity during the firing event. The amount of fluid entrained by the barb made a noticeable difference in the speed of the barb after the initial acceleration phase.

Our results indicate that although the initial acceleration phase is sufficient to transition the barb into an inertial regime allowing the barb to reach the target, it is likely other components of the firing mechanism are necessary to sustain the inertial environment for the barb to complete its course. In our previous model where the barb moves at a prescribed velocity, the effects from the boundary layer interactions only affect the prey rather than truly capturing an interaction. Results from [21] showed that the size of the prey was most important in resisting the effects of the approaching fluid front. In the current model, the length of the barb becomes more important as the two structures interact, with a longer barb able to sustain higher velocities for longer periods of time. For the longer barb, the boundary layer is narrower and longer, and this enhances the barb’s ability to reach the prey by accelerating the surrounding fluid and creating a layer that acts as a momentum reservoir. At the same time, the boundary layer also acts as a impedance as the barb nears the prey. As the boundary layer from the prey interacts with the moving fluid front of the barb, the barb is slowed while the prey is pushed out of the way.

Beyond the high accelerations needed to reach a Reynolds number sufficient for the barb to successfully reach the prey, high initial accelerations are needed for the barb to reach a velocity upon impact sufficient to puncture the cellular membrane or other external covering of the prey. This velocity is important since the likelihood of puncture is proportional to the square of barb’s velocity as the kinetic energy of the barb is transferred to the prey upon contact [36]. High initial accelerations alone, however, may not sustain the velocity of the barb for long enough to puncture a prey. These results suggest that other mechanisms, such as high osmotic pressures gradients, could be required for the overall mechanism to be successful. It is known that high osmotic pressures are involved in nematocyst firing [37]. The mechanism itself involves multiple stages including the initial release and the subsequent injection of the tubule. High osmotic pressures could create a sustained, accelerating front of fluid that may further serve to push the barb towards the prey.

In the future, we will investigate the boundary layer interactions in more detail, including the effects of relative motion, commonly observed strategies such as multiple nematocysts firing, and the presence of a background flow velocity due to osmotic pressure. Ultimately we plan to move our studies to three dimensions to more fully examine the mechanism involved in small-scale ballistics across Reynolds numbers. We have not yet incorporated puncturing dynamics into this model; we instead measure the final velocity of the barb at the time of contact with the prey (if it occurs) as an indication of whether puncturing is possible. Puncture dynamics represent a significant separate modeling challenge, and will be part of future work.

Finally, a natural extension of this work is to move into three-dimensions. We note that the current two-dimensional simulations represent an accelerating plate moving towards a cylinder of infinite length in the third dimension. Furthermore, vortex-stretching which may be relevant at higher $Re$ is not captured in two dimensional simulations. We anticipate that the same general trends will be observed in three dimensions, namely that the barb experiences rapid deceleration and that prey may be pushed out of the way at low $Re$. The 2D simulations can be used to develop a targeted 3D study that is focused on $Re$ between 1 and 100 and the length of the barb.