Quantification of the influence of drugs on zebrafish larvae swimming kinematics and energetics

Ethics

Animal work was carried out in compliance with the Animal Ethics and Welfare Committee, Department of Life Sciences, Glasgow Caledonian University, and UK Home Office under Project License PPL 60/4169.

Experiment setup

The experiment setup (Fig. 1A) and main methodologies were developed in our previous study for the determination of the toxicity of acrylamide on zebrafish locomotion via a colour preference experiment (Jia et al., 2017). In the present study, 85 dpf (days post-fertilization) wildtype zebrafish (Danio rerio) siblings were divided into four groups: the first group (control group) including 20 larvae was immersed in E3 medium (five mM NaCl, 0.17 mM KCl, 0.33 mM CaCl₂, 0.22 mM MgSO₄, and 0.1% methylene blue); the second group including 20 larvae was immersed in E3 medium with 0.01% acetic acid for 10 min. The third group including 20 larvae was immersed in 500 μM DPH solution (5, 5-DPH sodium salt) and the fourth group including 20 larvae was immersed in 100 mg/L yohimbine solution (yohimbine hydrochloride). The petri dish was illuminated by a light-emitting diode panel, driven by an adjustable DC power supply (CSI5003XE; Circuit Specialists, Tempe, AZ, USA) to provide a continuous and constant light. A high-speed video camera (EoSens CL MC1362; Mikrotron, Unterschleißheim, Germany) was used to record fish swimming behaviour. The frame rate of the camera was set at 500 frames per second during the entire experiment process. As in the subsequent CFD numerical modelling, the selected fish with a tail beating frequency being less than 70 Hz, thus, 7–8 frames within one beat cycle is sufficient to capture the fish tail motion. The water temperature was set at as 27 °C as this is the common temperature widely used in zebrafish experiments. Before the camera started recording, fish in all groups were allowed to swim freely for about 10 min to adapt to the water environment, once recording started, there was no stimulation to force the fish swimming forward. In this study, only quasi-steady cruising swimming regime is investigated excluding the sudden-start process. This is evident from previous research that cruising with cyclically motion is essential for fish larvae to cover the distance for migration and dispersal (Sancho, Ma & Lobel, 1997). In addition, cruising swimming has been studied extensively, which makes it easier to compare with other researchers’ results.

Data processing algorithm

An in-house MATLAB code was developed and used to post-process the recorded videos and extract zebrafish swimming kinematic characteristics, that is motion equations. Figure 1B depicts the key steps for the process. The original image recorded from the camera was converted to a binary image consisting of the sketch of zebrafish larva only with ‘im2bw’ function in MATLAB image processing toolbox. With some adjustments and ‘bwboundaries’ function in MATLAB, a binary image of zebrafish can be extracted, the entire position vector can be obtained for points distributed on fish outline. All images were skeletonised into a single backbone curve using functions ‘bwmorph’ and ‘thin’ operation.

The coordinated pixels on the backbone curve were then divided into equal-distant segments. These segments were simplified as connected straight lines to calculate relative orientation variation with time between two adjacent segments using MATLAB curve fitting toolbox. Physical representation of the intersection angle is shown in Fig. 1C and the mathematical intersection angle is expressed by Fig. 1D, and calculated with Eq. (1), where i denotes the points numbering from one, and θ_j is the relative angle between each two body segments. As elucidated by Muller (2004), the travelling wave of curvature travels along the fish body at a near constant rate, thus an averaged frequency was selected for the entire relative orientation functions. Equation (2) represents a sample prescribed motion equation for relative angle θ_j between each two body segments.

CFD solver and motion solver

Zebrafish larva CFD model

Zebrafish larva model used in OpenFOAM (22) (https://www.openfoam.com/) was built with 51 ellipses extracted from the real fish silhouette and controlled by nine deformation equations as shown in Fig. 2A. To simplify the model, the eyes and fin fold are excluded in the CFD fish model. Density of the fish is assumed to be the same as water. The same average body parameters such as body segments mass and length listed in Table 1 is used for all fish in CFD simulation. The flow field was numerically simulated using the open source CFD toolbox OpenFOAM version 3.0.x. The 3-D computational domain is 15 times the fish body length in the longitudinal (x) direction, 10 times of fish body length in transverse (y) direction and four times of fish body length in perpendicular (z) direction. The overall fluid domain is assumed to be at rest initially. In the simulation, the medium is water, therefore, the kinematic viscosity of the fluid υ, which can be expressed as ${\displaystyle \frac{\mathrm{\mu }}{\mathrm{\rho }}}$, is 10⁻⁶ m²/s. Pressure boundary conditions are taken as zero gradient for all boundaries except the front and back plane, which are set as symmetry; velocity boundary conditions for fish model were taken as moving wall velocity for all body segments and fixed value for the remaining patches. Local mesh around fish model were depicted in Fig. 2B. Mesh around head region before the vertical dash line shown in Fig. 2C is enlarged to be clearer. Considering the constraints in OpenFOAM regarding large mesh deformation to model self-propelled zebrafish swimming, fully unstructured mesh was used to tolerate the internal mesh deformation. For ellipses with high aspect ratio, the mesh is specially refined at the tips to ensure that enough cells are present to precisely capture the vortex at the tips. The Reynolds number is defined as ${\displaystyle \frac{vL}{\mathrm{\upsilon }}}$, v stands for final constant forward swimming velocity, L is the body length of fish larva, and υ represents the kinematic viscosity. For the entire simulation, Reynolds number is set as 300, which stands for intermediate flow regime.

Table 1:

Detailed zebrafish larva body information.

Body section number	Mass (mg)	Length (mm)
1	0.0385	0.4
2	0.0553	0.4
3	0.0425	0.48
4	0.0308	0.48
5	0.0212	0.48
6	0.019	0.48
7	0.011	0.48
8	0.009	0.48
9	0.003	0.32
Total	0.23	4

DOI: 10.7717/peerj.8374/table-1

Validations for multi-body coupling

Grid independence and sensitivity test

A grid independence test was carried out on a self-propelled 0.01% acetic acid treated zebrafish model with fully prescribed deformation with two mesh sizes, medium and fine mesh. In addition, the time step influence is also tested, for example Case 1: Medium mesh (M) with 83,200 total cells and time step of T/650, Case 2: Medium mesh and smaller time step (MS) of 83,200 cells and time step of T/1300, and Case 3: Fine mesh (F) of 166,400 cells and time step of T/650. The size of the fluid domain remains same for three cases.

To exclude mesh resolution and time step size influences on simulation results, we have compared the forward velocity and total force of zebrafish larvae treated with 0.01% acetic acid and the results are shown in Figs. 3A and 3B. Computational results for three cases including kinematic performance and fluid domain calculation shows close results. To save computational time and keep accuracy, we use the mesh formation and time step size the same as Case 1.

In addition, we have also tested the sensitivity of the results to body segments of zebrafish larvae model with CFD simulation. Considering the accuracy and efficiency for body segmentation, we have divided the body trunk into 5, 10 and 15 segments, respectively. By capturing the body deformation and simulate the forward motion with CFD toolbox, we compared the forward velocity and total hydrodynamic force in Figs. 3C and 3D, indicating that the simulation results are not sensitive to the number of body segments.

Multi-body dynamics validation

To validate our numerical coupling methodology on multi-body structure simulation, a comparison based on jellyfish-inspired swimming provided by Wilson was performed (Wilson & Eldredge, 2011). Figure 4A depicts the jellyfish shape, points on the figure indicate how the main body is separated. Figures 4B and 4C illustrate the mathematical model we have created and our CFD model. The prescribed angle functions are represented by θ₁, θ₂ and θ₃. Detailed functions are expressed with Eq. (6). By specifying relative orientation equations between every two sections, the jellyfish model can move upwards with alternant contraction and refilling. Different Reynolds numbers were tested, which is defined by the maximum jellyfish diameter and undulation period, that is ${\displaystyle \frac{{D^2}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{T\mathrm{\upsilon }}}$, D_max is the maximum diameter of the jellyfish model, T is the undulation period, and υ denotes kinematic viscosity. Here we selected Re equivalent to 140 and 70. Figures 5A and 5B depict the kinematic performance of the multi-body structure with fully prescribed motion, while Fig. 5C illustrates the power required for the model to move cyclically and Fig. 5D shows the vorticity comparison at Re = 140. All of our results are in close consistent with Wilson & Eldredge’s results, the small discrepancies can be caused by the different numerical methods used.

(6) $$\begin{array}{l}{\mathrm{\theta }}_1=-0.1472\cos \left(6.538t\right)-0.3247\sin \left(6.538t\right)+0.7551\\ {\mathrm{\theta }}_2=-02366\cos \left(6.456t\right)+0.1645\sin \left(6.456t\right)+0.3472\\ {\mathrm{\theta }}_3=-0.08218\cos \left(6.427t\right)+0.04095\sin \left(6.427t\right)+0.4511\end{array}$$

Validations on numerical and experimental methodologies

The first objective of our study was to verify the correctness/accuracy and consistence of proposed experiment and CFD method in this work. To achieve this goal, we used the time-varied angle data, extracted via post-processing of experiment video records as input data to CFD fish body discrete elements, the forward motion of fish is the solution by solving a coupled flow solver and multi-body dynamic method as described in “CFD Solver and Motion Solver”. The first and last fish body segment motion is not prescribed whereas they are our numerical solutions, which represent the fish head and tail locomotion. With simulation results, we first compared our yawing head angle and tail beat angle with data collected from observed zebrafish swimming. Figures 6A and 6B show an average head and tail angle for 10 individual zebrafish larvae. As seen from the figures, both the head and tail-beat angles were calculated based on the prescribed constrained deformation equations of zebrafish-matched experimentally observed data. Slight differences can be caused by subtle insufficient accuracy of the fish model capturing toolbox, which might lead to slight inconsistency in certain captured points of head/tail angle. The forward swimming velocity shown in Fig. 6C reaches about 95% of the experimentally measured swimming speed, approximately 19.25 body lengths per second. The small gap exists as the real fish may not bend the whole body as symmetrically as our model, which is fully controlled by prescribed sinusoidal functions. The subtle modelling errors may additionally lead to differences in swimming speed. We have also calculated the Strouhal number (St) for simulated CFD results. The St is defined as ${\displaystyle \frac{Af}{U}}$, where A is the tail beat peak-to-peak amplitude, f is the tail beat frequency and U is the averaged swimming velocity. In our simulation, the St is approximately 0.8. A statistical data about forward velocity for 10 zebrafish larvae is shown in Fig. 6D (to save space of the image, only ten results are shown).

Figure 7 depicts the vorticity iso-surfaces formed based on Q Criterion behind a swimming normal zebrafish larva at different instants in time within one time period and the dorsal view for vorticity iso-surfaces. Here, Q can describe the wake topology and defines vortices as positive second invariant of velocity gradient in region where vorticity magnitude is greater than strain-rate magnitude (Kolář, 2007). As seen from the most left and right column of Fig. 7, flow patterns behind the fish are represented by detached vortices and shown as translucent green fragments. Vortices starts to form in the vicinity of head, transmits downstream to tail and detaches at the tail, which are consistent with the fish tail motion; when the lateral displacement of the tail reaches the highest amplitude, vortices starts to shed at the tail tip, the already formed vorticity in the wake are mixed with the newly formed vorticity at tail tip. The most right column also displays a 3-D view of the vortex rings generated behind the fish to understand formation of flow patterns better. To validate the numerical methodology, Fig. 7 also compares the body curvatures of CFD model and the real fish in the recorded experiment video. As can be seen, two sets of results match very well in terms of body shape at all specific time within a period, indicating that our CFD model is able to imitate the self-propelled swimming of zebrafish larva and its interactions with surrounding fluid.

Power distribution along fish body and an initial approximation of the cost of transport

As the movement of each two neighbouring body segments is constrained with a prescribed deformation equation except for fish head and tail, mechanical power distribution along the fish body can be approximated by power generated by each body segment. The mechanical power generated from fish muscle includes the translational power due to linear motion and the rotational power due to body rotation. As the fish is moving cyclically, all the other terms are cancelled out except for the rotational power. Therefore, the mechanical power is estimated with the cross product of torque and angular velocity shown as Eq. (7a), and the total power transmitted into the water is calculated with Eq. (7b).

In the above equation, $P_M$ is the mechanical power of fish muscle, and $P_H$ represents hydrodynamic power generated by interactions with surrounding fluid. M_i is the internal torque for the ith joint calculated by MBDyn in the global frame, ω_i represents the angular velocity for the ith joint. F_j is the hydrodynamic force acting on the jth body and V_j represents the jth body velocity.

During muscle contraction, fish body bends, energy is generated and transmitted into the water, and the bended body interacts with the surrounding fluid, a thrust force generates and pushes the fish moving forward. Within this process, Approximately 20% of the energy (depending on fish species, the percentage can float dramatically) (Zhang, Yu & Tong, 2014) generated by muscle contraction is consumed due to the viscous dissipation of fish body tissues, and the remaining energy is transmitted into the water. In our model, viscous dissipations of fish body tissues are neglected to simplify the simulation. Therefore, as the variation of kinetic energy is zero during cyclic swimming, the mechanical power is fully transformed into the hydrodynamic power, which means the absolute value hydrodynamic power equals the mechanical power, and we only need to show hydrodynamic power for 20 zebrafish larvae as depicted in Fig. 8B. The calculated hydrodynamic power has been changed into absolute values as the sign is different from mechanical power. To further understand the different power generated at different locations along the fish body, we have compared the time-history of forces and velocity in Fig. 8 at three typical points shown in Fig. 8A, representing head region, body region and tail region. The trajectory is not completely parallel to X axis in global frame, force and velocity shown in Figs. 8C and 8D are pointing towards the real moving direction of zebrafish. An approximation of hydrodynamic power distribution on each body section is also depicted. As shown in Fig. 8E, the averaged hydrodynamic power for 20 fish larvae shows a significant higher value starting from approximately 75% of body length. According to motion equations, this region has the largest motion amplitude along the body in global frame, resulting in larger fluid force, thus more hydrodynamic power. In Fig. 8F, the mechanical power generated along the body shows an increase towards the tail and a steep decrease at the tail.

Detailed mechanical and hydrodynamic power and an approximation of cost of transport are summarised in Table 2. The cost of transport is defined as energy spent to travel unit distance per unit mass, which is expressed as ${\displaystyle \frac{P_m}{U}}$, P_m is the power per unit mass, and U is the forward velocity. Limited by the size of the table, only 10 fish data was selected and displayed in the table.

Comparison of kinematics and energetics between normal group and drug treated groups

Swimming performance of zebrafish larvae might be affected by drugs with different concentrations. According to Lopez-Luna’s research (Lopez-Luna et al., 2017b), zebrafish larvae exposed to 0.01% acetic acid displayed more active responses than the normal zebrafish, and these active behaviours were sustained for longer. Liu et al. (2016) has tested influences of (DPH) on 5 dpf zebrafish larvae locomotion at different concentrations and found that exposure to higher concentrations of DPH under light condition leads to decreased locomotor activities. Their team also tested the effect of 100 mg/L yohimbine on 5 dpf zebrafish larvae, but no obvious effect has been found on locomotor activity under light condition (Li et al., 2015). By simulating the swimming behaviours of zebrafish larvae under these circumstance, we examined the forward swimming velocity and cost of transport differences after exposing to 0.01% acetic acid, 500 μM DPH, and 100 mg/L yohimbine. Figures 9A and 9D depict the averaged head and tail comparison for three drugs and control group zebrafish larvae. As the different frequencies and initial angles were used, phase differences existed among those groups in head and tail angle respectively. From the figure, there are no significant differences for the amplitude of head angle, whereas the maximum tail beat angle for zebrafish larva immersed in the 0.01% acetic acid showed a larger value compared with the rest groups. Figure 9B illustrates the forward swimming speed comparison between control group zebrafish and drug treated groups. The forward velocity for acetic acid treated zebrafish is dramatically higher than other groups’ zebrafish, but there is no significant differences between control group and yohimbine treated group. For zebrafish treated with 500 μM DPH solution, there is a significant decrease of velocity compared with other groups. These results indicate that the forward swimming velocity might be influenced by the body undulation frequency and tail beat amplitude. With larger frequency and tail beat amplitude being represented by maximum tail beat angle, the forward swimming speed tended to increase.

As the tail beat frequency and amplitude are both increased, energy generated by fish body will increase as well, which might give rise to lower efficiency as the side oscillations consumes more energy without contributions to thrust. Figure 9C displays comparisons of the hydrodynamic power P_H, cost of transport were calculated for all cases, resulting in 81.73 μJ/m·kg, 96.24 μJ/m·kg, 65.44 μJ/m·kg and 82.32 μJ/m·kg for control, acetic acid treated, DPH treated, and yohimbine treated group, respectively. These values are similar to those reported by Li et al. (2016) (from 105 μJ/m·kg to 50 μJ/m·kg) in on larval zebrafish. It appears that the DPH treated group performed with greater efficiency than the other groups as it has the lowest tail beat amplitude and frequency. Figure 9E depicts the forward velocity for all groups. Based on the assumption of no differences among those groups, we have calculated the P value, which is usually evaluated in statistical hypothesis testing to determine the reliability of the results. The P value is smaller than 0.0001 for control group comparison with DPH treated group and acetic acid treated group, which indicates that significant increment of forward velocity can be seen after 0.01% acetic acid treatment, and a decrease of velocity happens after treated with 500 μM DPH solution. There is no obvious differences between control group and yohimbine treated group, indicating that 100 mg/L concentration of yohimbine will not influence the locomotion behaviour dramatically of 5 dpf zebrafish larvae. Figure 9F displays similar comparison as Fig. 9E, but with the parameter changed to hydrodynamic power. The resulting P value is smaller than 0.0001 for control group comparison with DPH treated group and acetic acid treated group, leading to the conclusion that averaged hydrodynamic power of acetic acid treated zebrafish larvae is higher than the control group specimens, whereas the zebrafish larvae treated with DPH solution resulted in lower averaged hydrodynamic power. There is no obvious difference between control group and yohimbine treated group as well, which is consistent with results shown in Fig. 9E. Figure 10 compares the vorticity on x-y plane of the 3-D fish model. By comparing them in one period of time, it can be seen that the vortex detached from the tail tip is faster for acid treated group compared with control group fish. The earlier detached vortex has been labelled with black circles in the right column, indicating the larger distance travelled within one period of time, that is higher velocity. For DPH treated group, within same one period of time, the evolution of vorticity is much shorter than other groups, approximately half of other groups, suggesting that only half of the distance travelled by DPH treated group compared with control group. For yohimbine treated group, vorticity patterns are similar to those in control group. All of the vorticity results for drug treated groups are in consistent with velocity comparisons depicted in Fig. 9.